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ABSTRACT 


A  reduced  shallow  water  model  under  constant,  non-zero  advection  in  infinite  do¬ 
mains  is  considered.  High-Order  Givoli-Neta  (G-N)  and  Hagstrom-Hariharan  (H-H)  non¬ 
reflecting  boundary  conditions  (NRBCs)  are  introduced  to  create  a  finite  computational 
space  and  solved  using  a  spectral  element  formulation  with  high-order  time  integration. 
Numerical  examples  are  used  to  demonstrate  the  synergy  of  using  high-order  spatial,  time 
and  boundary  discretizations.  Several  alternatives  are  also  presented  for  solving  open  do¬ 
main  problems.  These  alternatives  include  adjustments  to  the  G-N  NRBC  based  on  phys¬ 
ical  arguments  as  well  as  formulating  the  boundary  condition  for  arbitrary  domains  using 
unstructured  grids.  The  H-H  polar  NRBC  is  also  formulated  in  an  unstructured  grid  setting 
and  extended  to  include  dispersive  effects.  Results  show  that  by  balancing  all  numerical 
errors  involved,  high-order  accuracy  can  be  achieved  for  unbounded  channel  problems. 
Further,  the  adjustments  to  the  G-N  and  H-H  NRBCs  to  operate  in  an  unstructured  grid 
setting  are  shown  to  significantly  reduce  errors  over  first  order  non-reflecting  boundary 
schemes  when  operating  in  an  open  domain  configuration. 
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I.  INTRODUCTION 


Simulation  of  wave  propagation  through  large  -  perhaps  unbounded  -  domains  has 
been  an  active  area  of  research  for  several  decades.  Such  studies  are  important  to  many 
applications  such  as  acoustics,  electromagnetics,  meteorology,  solid  geophysics  and  aero¬ 
dynamics  to  name  just  a  few.  The  combination  of  the  complexity  of  the  partial  differential 
equation  sets  involved  and  the  infinite  possibilities  of  initial  data  mandates  the  numerical 
solution  to  such  problems.  Of  course,  to  undertake  the  numerical  solution  on  an  infinite 
domain  would  be  both  foolhardy  and  impossible.  Generally  speaking,  to  overcome  this 
computational  challenge,  it  is  quite  common  to  truncate  the  infinite  domain  by  imposing 
some  type  of  boundary  condition  on  a  “sufficiently  large”  sub-domain  that  captures  the  area 
of  interest. 

When  truncating  the  domain,  the  modeler  must  devise  boundary  conditions  for  the 
truncated  domain.  Of  course,  by  imposing  a  boundary  where  one  does  not  physically 
exist,  the  problem  is  changed  -  and  unless  chosen  carefully,  would  certainly  be  expected  to 
pollute  the  solution  as  the  problem  evolves  and  impinges  on  the  boundary.  Therefore,  two 
main  possibilities  exist  for  the  modeler: 

•  Choose  a  convenient,  easily  implementable  boundary  condition  that  does  not  nec¬ 
essarily  reflect  the  physical  problem  and  solve  it  on  a  large  sub-domain.  The  idea 
behind  this  technique  is  that  the  boundary  effects  are  negligible  for  a  short  time  evo¬ 
lution  of  the  problem  in  a  small  area  of  interest  away  from  the  boundaries. 

•  Choose  a  boundary  condition  that  preserves  the  true  behavior  of  the  infinite  solution 
at  the  boundary  and  solve  the  problem  on  a  smaller  sub-domain.  The  idea  behind  this 
technique  is  that  the  additional  effort  extended  to  impose  a  better  boundary  condition 
will  be  worth  the  effort  and  allow  for  solving  the  problem  on  a  smaller  domain. 

For  obvious  reasons,  the  first  possibility  has  only  limited  usefulness.  To  see  why,  suppose 
that  we  wanted  to  model  the  wave  motion  following  a  pebble  dropped  in  the  center  of  a 
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large,  still  pond.  Now,  suppose  that  we  have  a  truncated  domain  to  model  this  phenomena 
-  say  a  bathtub.  If  the  pebble  is  dropped  in  the  bathtub,  the  waves  generated  by  the  pebble 
would  propagate  much  like  that  in  the  pond  -  until,  that  is  -  the  wave  front  reaches  the  hard 
walls  of  the  bathtub.  At  this  point,  the  bathtub  model  ceases  to  be  a  useful  representation 
of  the  pond  due  to  behavior  caused  by  the  non-physical  boundary.  If  the  modeler  wishes  to 
see  what  happens  a  short  time  later  -  a  larger  bathtub  would  be  required.  This  same  prin¬ 
ciple  would  apply  for  the  numerical  solution  of  this  propagation  problem  -  a  poor  choice 
of  boundary  condition  mandates  the  use  of  a  larger  computational  domain.  This  in  turn, 
requires  additional  computational  resources.  For  this  reason,  much  effort  has  and  continues 
to  be  exerted  on  finding  suitable  boundary  conditions  that  apply  on  smaller  domains. 

This  dissertation  examines  the  use  of  high-order  non-reflecting  boundary  conditions 
(NRBCs)  to  solve  a  class  of  infinite  domain,  wave  propagation  problems.  In  the  last  35 
years  or  so,  much  research  has  been  done  to  develop  NRBCs  that,  after  discretization,  lead 
to  a  scheme  that  is  stable,  accurate,  efficient  and  easy  to  implement.  Of  course,  it  is  difficult 
to  find  a  single  NRBC  that  is  ideal  in  all  respects  and  all  cases;  this  is  why  the  quest  for 
better  NRBCs  and  their  associated  discretization  schemes  continues. 

Sequences  of  increasing-order  NRBCs  have  been  available  for  a  long  time  (e.g.,  the 
Bayliss-Turkel  conditions  [1]  constitute  such  a  sequence),  but  they  had  been  regarded  as 
impractical  beyond  2nd  or  3rd  order  from  the  implementation  point  of  view.  Only  since  the 
mid  90s  have  practical  high-order  NRBCs  been  devised.  The  first  high-order  local  NRBC 
was  proposed  by  Collino  [2],  for  two-dimensional  time -dependent  waves  in  rectangular 
domains.  Its  construction  requires  the  solution  of  the  one-dimensional  wave  equation  on 
the  boundary.  Grote  and  Keller  [3]  developed  a  high-order  converging  NRBC  for  the  three- 
dimensional  time -dependent  wave  equation,  based  on  spherical  harmonic  transformations. 
Sofronov  [4,  5]  proposed  exact  boundary  conditions  for  the  three-  and  two-  dimensional 
wave  equations  in  spherical  and  polar  coordinates,  respectively  (it  is  proved  that  NRBCs 
demonstrated  in  [3]  and  [4]  are  reduced  to  each  other). 
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Hagstrom  and  Hariharan  [6]  constructed  high-order  NRBCs  for  the  two-  and  three- 
dimensional  time-dependent  wave  equations  based  on  the  analytic  series  representation  for 
the  outgoing  solutions  of  these  equations.  For  time-dependent  waves  in  a  two-dimensional 
wave  guide,  Guddati  and  Tassoulas  [7]  devised  a  high-order  NRBC  by  using  rational  ap¬ 
proximations  and  recursive  continued  fractions.  Givoli  [8]  has  shown  how  to  derive  high- 
order  NRBCs  for  a  general  class  of  wave  problems,  leading  to  a  symmetric  FE  formulation. 
In  [9],  this  methodology  was  applied  to  the  particular  case  of  time-harmonic  waves,  using 
optimally  localized  Dirichlet-to-Neumann  (DtN)  NRBCs  (see  also  [10]). 

Previous  studies  of  this  nature  have  encountered  limits  in  the  accuracy  of  such  solu¬ 
tions.  These  accuracy  limits  can  be  caused  by  time  and  space  discretization  as  well  as  from 
the  boundary  scheme  used  in  the  solution  of  the  problem.  We  seek  to  employ  high-order 
numerical  methods  in  time  and  space  to  diminish  the  effects  of  discretization  error  in  order 
to  determine  the  true  efficacy  of  a  given  non-reflecting  boundary  condition. 

The  rest  of  the  dissertation  is  outlined  as  follows.  Chapter  II  motivates  and  derives 
the  equations  for  the  problem  under  consideration.  In  Chapter  III,  we  summarize  the  main 
boundary  schemes  currently  in  use  and  specifically  show  the  Givoli-Neta  (G-N)  NRBC 
in  detail.  We  then  describe  the  high-order  spectral  element  method  used  to  discretize  the 
problem  in  space  (up  to  16th  order  polynomials)  in  Chapter  IV.  Chapter  V  discusses  the 
Runge-Kutta  time  discretization  demonstrated  up  to  lQth  order.  Chapter  VI  provides  nu¬ 
merical  examples  in  various  configurations  and  conditions  to  demonstrate  concepts.  Chap¬ 
ter  VII  considers  challenges  associated  with  arbitrary  boundary  configurations  and  provides 
results  for  a  low-order  boundary  treatment  using  unstructured  grids.  Chapter  VIII  exam¬ 
ines  the  potential  for  exploiting  other  non-reflecting  boundary  conditions  employed  in  an 
unstructured  grid  formulation  and  concludes  with  areas  of  future  research. 
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II.  EQUATIONS  OF  FLUID  MOTION 


In  order  to  scope  the  enormous  problem  of  wave  propagation  and  provide  a  good 
test  bed  of  examples  for  simulation  of  the  non-reflecting  boundary  conditions  under  exami¬ 
nation,  we  will  derive  the  shallow  water  equations.  The  set  of  equations  under  consideration 
have  been  used  to  predict  Tsunamis  and  storm  surges  [1 1],  as  well  as  modeling  atmospheric 
flows.  The  term  “shallow  water”  is  a  bit  deceiving,  as  the  medium  does  not  necessarily  have 
to  be  water,  nor  does  it  have  to  be  shallow.  To  elaborate  -  the  equations  under  consideration 
are  generated  by  very  general  physical  principles,  namely  conservation  of  mass  and  mo¬ 
mentum,  that  are  then  simplified  using  reasonable  assumptions.  In  this  derivation,  we  will 
abbreviate  the  approach  taken  by  Cushman-Roisin  [12],  Pedlosky  [13]  and  Batchelor  [14] 
and  further  simplify  the  equations  by  reducing  them  to  a  scalar  Klein-Gordon  equation 
equivalent. 

A.  CONSERVATION  OF  MASS 

One  of  the  fundamental  physical  principles  is  that  mass  can  neither  be  created  nor 
destroyed.  Consider  a  control  volume  -  a  fixed  region  in  space  where  fluid  is  allowed 
to  occupy  and  pass  through.  Within  this  control  volume,  mass  is  conserved.  In  other 
words,  for  the  mass  to  change  in  a  control  volume,  there  must  be  mass  passing  through  the 
boundary  of  the  control  volume.  Suppose  dm  is  an  infinitesimal  portion  of  the  mass,  dV 
and  dA  are  infinitesimal  portions  of  the  control  volume  and  its  boundary,  respectively,  and 
p  is  the  density  of  the  fluid  occupying  the  control  volume.  Then  by  this  argument,  the  mass 
enclosed  by  the  surface  at  any  instant  is  f  dm  =  f  p  dV .  Further,  the  net  rate  that  mass 
is  flowing  outwards  across  the  boundary  is  f  (p u)  •  n  dA.  Putting  these  two  ideas  together 
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we  find  that  for  mass  to  be  conserved: 


I  pdV  =  -  J  (pu)  •  n  dA. 

Time  rate  of  mass  Net  rate  of  mass  flux 

change  in  CV  across  boundary 

Here:  u  =  (u.  v.  w)T  is  the  fluid  velocity  and  n  is  the  outward  pointing  unit  normal  on  the 
boundary. 

Further,  upon  differentiation  under  the  integral  sign  (remembering  that  the  control 
volume  is  fixed  in  space)  and  transforming  the  surface  integral  using  the  divergence  theo¬ 
rem: 


/  ^  dV  +  J  V  ■  (pu)  dV  =  0 
/(|+v.pu)  dV=0 

This  relation  is  valid  for  all  choices  of  the  control  volume,  therefore,  the  integrand 
must  also  be  identically  zero,  i.e., 

^  +  V  •  pu  =  0  (III) 

The  differential  equation  (II.  1)  is  commonly  referred  to  as  the  continuity  equation  in  fluid 
mechanics. 

B.  IMPORTANCE  OF  THE  EARTH’S  ROTATION 

Modeling  phenomena  on  a  large  scale,  such  as  the  currents  of  the  ocean  or  winds  in 
the  atmosphere,  may  require  special  handling  since  the  earth  is  not  static.  In  fact,  the  earth 
rotates  on  its  axis  approximately  once  every  24  hours.  This  results  in  a  mean  rotation  rate 
fl  of: 

n  = _ ^ _ . 

rotation  period 
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Figure  1:  Fixed  (A",  Y)  and  rotating  (x,  y )  frameworks  of  reference. 

When  the  “exact”1  rotation  period  of  the  earth  is  considered,  this  results  in  a  mean  angular 
velocity  of  7.292  x  10~5ra^ns  [15].  The  trajectory  of  a  fluid  in  motion  is  expected  to  be 
influenced  by  this  rotation  if  the  fluid  traveling  at  the  speed  U  covers  a  distance  L  in  a  time 
interval  greater  than  the  rotation  period.  This  concept  is  captured  by  a  non-dimensional 
parameter  £  and  is  defined  as: 

27t/0  27 tU 

£  =  L/U  =  ~VtL 

If  e  is  on  the  order  of  or  less  than  unity  (e  <  1),  then  we  would  expect  that  rotation 
is  important.  This  number  (neglecting  the  constant  multiple  2n)  is  known  as  the  Rossby 
number  [12]. 

1.  Equations  in  a  Rotating  Frame 

Now,  we  examine  this  rotating  frame  of  reference  on  the  earth.  From  the  human 
perspective  on  the  surface  of  the  earth,  we  appear  to  be  on  a  2-dimensional  surface.  Suppose 
that  we  have  X—  and  Y  —  axes  that  are  the  fixed  or  inertial  reference  frame  and  x —  and 
y—  axes  that  form  the  same  reference  frame,  but  rotating  at  the  angular  rate  of  Q.  The  unit 
vectors  that  follow  this  convention  are  defined  by  (I,J)  and  (i,j)  as  shown  in  Figure  l2.  It 

'The  mean  day  is  23  hours,  56  minutes  and  4.091  seconds,  but  variations  caused  by  friction  from  the 
earth’s  tides,  as  well  as  significant  geophysical  events  on  earth  have  been  observed  to  cause  fluctuations  in 
this  measure. 

2Adapted  from  [12],  Figure  2-1,  p.  17. 


7 


follows  that 


i  =  I  cos  fit  +  J  sin  f It  j  =  — I  sin  fit  +  J  cos  fit 

and  the  coordinates  of  the  position  vector  r  =  XI  +  Y  J  =  xi  +  yj  are  correspondingly 

x  =  X  cos  fit  +  y  sin  fit 
y  =  —X  sin  fit  +  Y  cos  fit 

Differentiating  once  with  respect  to  time  gives  the  rate  of  change  of  the  coordinates  relative 
to  the  moving  frame,  u  =  ^  =  ui  +  vj  (relative  velocity).  Differentiating  again  with 
respect  to  time  gives  the  rate  of  change  of  the  relative  velocity  in  the  moving  frame,  a  = 
^  =  ai  +  bj  (relative  acceleration).  When  completed  and  simplified,  the  absolute 

acceleration  in  the  inertial  frame  with  respect  to  the  relative  acceleration  is: 

A  =  (a  —  2f Iv  —  f i  +  (b  +  2flw  —  fl2r/)  j  (II. 2) 

=  Ai  +  Bj 

These  results  could  also  be  derived  in  a  vector  form  as  outlined  in  [13]  by  defining 
the  vector  rotation  in  a  direction  common  to  both  the  rotating  and  inertial  frames  of  refer¬ 
ence  -  f 2  =  flk,  where  k  is  a  unit  vector  orthogonal  to  the  plane.  In  this  case,  we  can  write 
(II. 2)  as: 

A  =  a  +  2flxu  +  flx(flxx)  (II.3) 

where  u  is  extended  to  u  =  wi+nj+tck.  It  should  be  noted  that  there  are  three  contributions 
to  the  acceleration  in  the  rotating  frame:  relative  acceleration  (a),  one  proportional  to  fl  and 
the  velocity,  and  one  proportional  to  fl2  and  the  position.  The  contribution  proportional  to 
fl  and  the  velocity  is  known  as  the  Coriolis  acceleration  and  the  other,  proportional  to  Q2 
and  the  position  is  the  centripetal  acceleration. 
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For  practical  purposes,  the  centripetal  acceleration  terms  are  often  neglected  since 
fl2  ~  O  (10-9).  Additionally,  even  though  centripetal  acceleration  causes  objects  on  the 
surface  of  the  planet  to  feel  an  outward  pull,  these  objects  do  not  fly  off  into  space.  In 
fact,  even  objects  at  rest,  (u,  v )  =  0  thus  removing  the  Coriolis  effect,  for  all  intents  and 
purposes,  remain  at  rest  as  the  gravitational  pull  of  the  earth  keeps  centripetal  acceleration 
in  check. 

When  we  neglect  the  centripetal  acceleration  terms  in  (II. 3),  the  absolute  accelera¬ 
tion  terms  in  the  inertial  frame  simplify  to 

A  =  a  +  2  fl  x  u 

=  (a  —  2fk>)  i  +  (b  +  20rt)  j  +  ck  (II.4) 

2.  3-Dimensional  Rotating  Earth  Model 

Now,  we  consolidate  our  results  to  apply  to  a  3-dimensional  rotating  earth  model. 
Consider  Figure  23  where  Q  is  oriented  along  the  axis  of  rotation  and  an  object  is  located 
on  the  surface  at  a  latitude  0.  A  local  coordinate  system  is  set  up  with  the  axis  orientation 
(x,  y,  z )  — *  (east, north, radial)  with  standard  convention  of  unit  normal  vectors.  In  this 
frame  of  reference,  the  earth’s  rotation  vector  is  expressed  as 

f2  =  cos  0j  +  fl  sin  0k.  (II. 5) 

Using  this  to  expand  (II.4),  we  find  that  the  acceleration  in  the  inertial  reference  in 
terms  of  the  rotating  components  has  the  following  components: 


i  : 

a  +  2 Qw  cos  0  —  2Qv  sin  0 

j  : 

b  +  2 Qu  sin  0 

(H.6) 

k  : 

c  —  2 Qu  cos  0. 

3  Adapted  from  [12],  Figure  2-8,  p.  27. 
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Figure  2:  Local  Cartesian  framework  on  spherical  earth. 

Here,  we  notice  the  terms  dependent  on  the  latitude  have  common  components,  namely: 

/  =  2f2  sin  0  (II.7) 

/*  =  2fi  cos  (f) 

The  coefficient  /  is  called  the  Coriolis  parameter  and  the  latter  is  the  reciprocal  Coriolis 
parameter.  If  we  examine  the  scales  of  the  components  described  in  (II. 6),  we  note  that 
if  we  are  describing  geophysical  fluid  motion,  the  depth  is  much  smaller  than  that  of  the 
surface  it  covers,  and  as  such,  the  flows  in  this  context  tend  to  be  parallel  to  the  surface 
minimizing  the  effects  of  vertical  flows.  For  our  model,  this  implies  that  the  effects  of  w 
are  negligible  compared  to  the  effects  of  u  and  v.  Additionally,  any  acceleration  induced  by 
the  rotation  in  the  vertical  direction  will  be  negligible  compared  to  those  along  the  surface. 
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This  simplifies  our  acceleration  in  the  inertial  reference  to: 


i  :  a  —  2Qv  sin  0 

j  :  b  +  2 £lu  sin  0 

k  :  c. 

In  vector  form,  this  can  be  written  as: 

A  =  a  +  /  (k  x  u)  (II. 8) 


Cushman-Rosin  [12]  provides  two  anecdotal  justifications  for  this  simplification.  While 
the  atmospheric  layer  that  determines  our  weather  is  only  about  10  km  thick,  cyclones  and 
anticyclones  spread  over  thousands  of  kilometers.  The  second  in  the  context  of  oceanic 
currents  notes  that  flows  are  generally  confined  to  the  upper  hundred  meters  but  spread 
over  tens  of  kilometers. 

C.  CONSERVATION  OF  MOMENTUM 


The  linear  momentum  of  an  object  of  mass  m  moving  with  velocity  u  is  defined 
to  be  the  product  of  the  mass  and  velocity:  P  =  mu,  and  in  a  closed  system,  must  be 
conserved.  Linear  momentum  is  related  to  the  forces  acting  using  Newton’s  second  law  of 
motion  -  namely,  that  the  time  rate  of  change  of  the  momentum  of  an  object  is  equal  to  the 
resultant  forces  on  the  object.  That  is, 


F  = 


(H-9) 


This  implies  that  if  resultant  forces  are  zero,  the  momentum  of  the  particle  must  be  constant. 
For  this  to  happen,  all  of  the  forces  (body  and  surface)  acting  on  an  object  must  sum  to  zero. 
In  terms  of  the  control  volume  framework  discussed  previously,  this  further  implies  that  if 
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momentum  is  entering  or  leaving  the  control  volume: 


d_ 

dt 


P 


d 

=  dt™* 


^  j  pwjW  =  j  pbdV  —  J  (pu)u-noL4  +  j  Tn dA. 


Body  forces 
acting  on  Q 


Net  momentum  flux 
across  bdry  of  fl 


Forces  acting 
on  bdry  of 


(II.  10) 


Here,  u*  is  the  velocity  in  the  inertial  frame.  Body  forces  are  those  acting  on  the  fluid 
volume  that  are  proportional  to  the  mass.  The  body  forces  considered  here  are  gravity  and 
(indirectly)  the  Coriolis  force  described  in  II. B.  Others  could  include  electromagnetic  and 
centrifugal  forces  pertinent  in  alternate  applications. 


1.  Gravity  Effects 

Gravity  acts  on  a  control  volume  strictly  towards  the  center  of  the  earth,  and  in  the 
local  coordinate  system  is  along  — k.  This  implies: 


pbg  =  -peg  k 


(11.11) 


Here,  g  is  the  gravitational  constant  which  varies  based  on  the  distance  from  the  center  of 
the  earth.  At  sea  level,  this  value  is  approximately  9.798(5,  and  in  this  context  can  be  taken 
to  be  constant. 


2.  Coriolis  Effects 

Coriolis  acts  in  a  way  to  “adjust”  the  control  volume’s  acceleration  when  going 
from  a  rotating  to  an  inertial  frame.  This  term  enters  via  the  expression  on  the  left  hand 
side  in  the  time  change  of  momentum.  Explicitly,  this  is 


Wt  f puJV  = 


dp  <9u* 

mu*  +  p^ 


dV=  I  (  (\i„  +PA)dV 


+  pa  +  pf  (k  x  u) 


dV 


d 


=  —  /  pudV  +  I  pf  (k  x  u)  dV 
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provided  that  p  is  constant  in  time. 


3.  Surface  Force  Effects 

Surface  forces  are  those  exerted  across  the  boundary  by  the  surrounding  matter. 
Typical  surface  forces  include  pressure  and  viscosity.  Since  we  are  examining  the  equations 
in  the  context  of  atmospheric  and  oceanic  applications,  the  effects  of  viscosity  are  small 
in  comparison  to  other  forces,  and  as  such,  will  not  be  considered  here.  VanJoolen  [16] 
derives  these  terms  in  detail  for  the  interested  reader,  but  further  concludes  the  contribution 
of  viscosity  can  be  neglected  in  this  application.  The  total  boundary  force  exerted  by  the 
surrounding  matter  at  the  surface  of  a  control  volume  is 

/  TncL4  =  —  [  pndA 


where  p  is  the  pressure  exerted  on  the  control  volume  by  the  surroundings.  This  surface  in¬ 
tegral  may  be  transformed  to  an  integral  over  the  volume  by  the  analogue  of  the  divergence 
theorem  for  a  scalar  quantity  [14,  p.  15]  giving 

—  J  pndA  =  —  j  VpdV. 


Alternate  derivations  of  this  quantity  can  be  found  in  [16]  and  [17],  and  provide  further 
insights  as  from  where  this  quantity  is  derived. 


4.  Momentum  Equations  in  a  Rotating  Frame 

We  now  apply  the  divergence  theorem  to  transform  the  surface  integral  in  (II.  10) 
and  collect  results. 


d_ 

dt 


pudV  + 


V  ■  (puu)  dV 


pg\tdV  + 


pf  (k  x  u)  dV  + 


VpdV  =  0 
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Here,  uu  is  a  tensor  product  of  the  velocity  vector.  Again,  differentiating  under  the  integral 
sign  and  consolidating  terms  yields: 


d  . 


V  •  (puu)  —  pgk  +  pf  (k  x  u)  +  Vp 


dV  =  0. 


Since  this  relation  is  valid  for  all  choices  of  the  control  volume,  the  integrand  must  identi¬ 
cally  be  zero,  i.e., 


—  (pu)  +  V  ■  (puu)  +  Vp  =  pgk  +  pf( u  x  k) . 


(11.12) 


Together  with  the  continuity  equation  (II.  1)  we  have  our  equations  of  fluid  motion  in  a 
rotating  frame. 

D.  FURTHER  SIMPLIFICATIONS  OF  THE  SYSTEM 


Our  equations  of  fluid  motion  have  taken  several  physical  principles  into  account, 
but  we  wish  to  simplify  them  more  in  order  to  get  a  tractable  model  that  can  serve  as  a 
launching  point  for  more  testing.  In  the  case  of  this  analysis,  we  wish  to  make  the  following 
assumptions  about  the  physical  problem  in  order  to  arrive  at  the  shallow  water  equations: 

The  fluid  is  homogeneous:  We  assume  that  the  fluid  density  is  constant  and  uniform  through¬ 
out  the  domain. 

The  fluid  is  inviscid:  This  implies  that  the  only  surface  force  acting  on  the  fluid  is  pressure 
(neglects  shear  forces  which  would  act  to  retard  the  motion  of  the  fluid.) 

The  fluid  is  incompressible:  Together  with  the  assumption  of  homogeneity  this  decou¬ 
ples  the  dynamics  from  any  thermodynamic  considerations  that  might  be  used  in 
another  setting. 

Centrifugal  forces  are  balanced  by  gravity:  This  allows  simplification  of  acceleration  terms 
in  the  rotating  frame  of  reference. 
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Figure  3:  The  shallow  water  model  with  irregular  bottom  topography. 

The  fluid  is  shallow:  The  depth  terms  in  the  applications  considered  are  much  smaller 
than  the  surface  it  covers,  and  therefore  implies  that  flows  are  primarily  along  the 
surface. 

We  consider  a  sheet  of  fluid  as  shown  in  Figure  3  with  properties  as  outlined  above. 
Here,  we  define  the  irregular  bottom  height  below  a  sensible  a  reference  value  z  —  0  as 
hs(x,y).  This  reference  level  could  be  considered  the  fluid  height  at  rest.  The  height  of 
the  surface  of  the  fluid  above  the  same  reference  level  is  defined  as  h(x,  y,  t).  The  depth 
of  the  fluid  is  therefore  H(x,  y,  t)  =  h(x,  y,  t )  +  hB(x,  y).  To  reiterate  the  shallow  water 
assumption,  we  have  a  value  for  the  height  of  the  fluid  II  (x.  y.  t )  that  is  much  smaller  than 
the  length  and  width  of  the  fluid.  Again,  we  consider  a  fluid  that  is  in  the  presence  of 
rotation  O  about  the  z— axis. 

1.  Mathematical  Simplifications 

As  outlined  in  [17,  Appendix  A],  the  momentum  equations  (11.12)  can  be  mathe¬ 
matically  simplified  using  nothing  but  product  rule  expansions  and  the  continuity  equation 
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to  yield: 


du  du  du  du 

ai  +  udx  +  vav 
dv  dv  dv 

ai  +  udx+vdy 

dw  dw  dw  dw 

m+uai  +  vav+wal 


W 


w 


dz 

dv 

dz 


1  dp 

pdx~  V 
1  dp 

pdy  =  ~fU 
1  dp__ 
p  dz  ^ 


(11.13) 


2.  Implication  of  Homogeneity 

Since  the  fluid  is  assumed  to  be  homogeneous  in  nature  (p  is  constant),  the  conti¬ 
nuity  equation  (II.  1)  reduces  to 


du  dv  dw 
dx  +  dy  dz 


V  ■  u  =  0. 


(II.14) 


3.  Implication  of  Shallow  Fluid 

This  assumption  allows  significant  simplification  of  our  fluid  motion  model.  Here, 
we  assume  that  the  surface  scale  of  the  problem  at  hand  is  much  larger  than  that  of  any 
depth  considerations.  Pedlosky  [13]  outlines  a  scaling  argument  that  shows  how  the  rel¬ 
ative  importance  of  terms  in  the  z— momentum  equation  allow  all  of  the  terms  except  for 
the  pressure  derivative  and  the  gravity  to  be  neglected.  This  collapses  the  z— momentum 
equation  to 


dp 

dz 


~ P9 ■ 


We  can  then  immediately  depth  integrate  to  yield 


V  =  -PQZ  +  A(x,y,t). 
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If  the  surface  is  under  some  constant  ambient  pressure  p(x,  y.  h )  =  p0,  this  implies  that 
A(x,  y,  t)  =  po  +  pgh(x,  y,  t)  thereby,  giving  us  an  expression  for  p 


P(x,  y,  h)  =  pg  (h(x,  y,  t)  -  z)  +  p0. 


We  compute  the  pressure  gradients  in  the  x—  and  y—  directions 


dp  dh(x,y,t)  dp  dh(x,y,t ) 

dx  ^  dx  dy  ^  dy 


(11.15) 


Here,  we  note  [13]  that  the  pressure  gradients  are  independent  of  z  so  that  the  horizontal 
accelerations  must  also  be  independent  of  z.  For  consistency,  we  therefore  assume  that  the 
horizontal  velocities  will  also  be  independent  of  z. 


a.  Primarily  Horizontal  Flow 

We  have  already  observed  that  since  the  flow  in  shallow  waters  is  primarily 
along  the  surface,  the  z—  momentum  collapsed  down  significantly.  We  can  use  this  argu¬ 
ment  to  similarly  simplify  the  x—  and  y—  momentum  equations.  In  this  case,  since  w  is 
significantly  smaller  than  u  or  v.  the  terms  and  w can  also  safely  be  neglected.  Sub¬ 
stituting  these  simplifications  along  with  our  pressure  gradients  (11.15)  into  (11.13)  results 
in  the  x—  and  y—  momentum  equations: 


du  du  du  dh(x,  y,  t ) 

m  +  udx  +  vdy  -  fv  =  -g^r- 
dv  dv  dv  dh(x,  y,  t ) 

m  +  udx  +  %  ~  fu  = 


(11.16) 
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b.  Continuity  Equation  Simplifications 

We  now  depth  integrate  our  continuity  equation  (11.14)  from  z  =  — hb(x ,  y  ) 
to  z  —  h(x ,  y,  t).  After  dropping  variable  dependencies  for  clarity,  this  yields: 


0 


V  ■  u  dz 


J  —hs 

f  (t  +  t) 

J-hB  \dx  dyJ 


dz  +  w \z=h  -  w\ z=-hg. 


(11.17) 


Considering  reasonable  boundary  conditions  for  the  last  term  as  detailed  in  [12],  we  specify 
no  normal  flow  on  the  rigid  bottom  (— hB )  and  a  corresponding  kinematic  condition  at  the 
fluid  surface  ( h ).  These  conditions  are: 


w(x,y,  -hB) 
w(x,y,  h,t) 


=  — u 

_  dh 
~  dt 


dhB 

dx 

+  u 


dhB 

V~dy~ 
dh  dh 
dx  dy 


(11.18) 


Combining  (11.17)  and  (11.18)  yields  (after  simplification  outlined  in  Appendix  A) 


dh 

dt 


d_ 

dx 


(Hu)  + 


d_ 

dy 


(Hu)  =  0. 


E.  LINEARIZING  THE  SHALLOW  WATER  MODEL 


(11.19) 


The  shallow  water  model  in  its  current  form  is  non-linear.  We  have  three  state  vari¬ 
ables,  u,  v,  and  h.  Each  of  these  are  defined  such  that  u  =  u(x,  y,  t)  and  v  =  v(x,  y,  t)  are 
the  unknown  fluid  velocities  in  the  x  and  y  directions,  h(x.  y,  t)  is  the  unknown  fluid  eleva¬ 
tion  above  the  reference  level,  /  is  the  Coriolis  parameter,  and  g  is  the  gravity  acceleration. 
We  now  introduce  the  following  shorthand  for  partial  derivatives 

„  _y_ 

^ CL  rx  ")  O db  rx  r\j 

oa  oaob 


18 


We  wish  to  find  a  linear  version  of  these  equations.  Right  now,  we  have 


dtu  +  udxu  +  vdyU  —  fv  =  —g  dxh 

dtv  +  udxv  +  vdyV  +  fu  —  —g  dyh  (11.20) 

d th  +  dx  (Hu)  +  dy  (Hv)  =  0. 

Now,  suppose  that  the  bottom  topography  is  flat  such  that  hB  is  constant  and  u  and  v  can 
be  described  by  a  constant  mean  term  and  a  small  0(5)  deviation  from  that  value,  i.e., 

u  =  U  +  u*  v  =  V  +  v*  H  =  Iib  +  h 

To  be  clear,  U  and  V  are  the  mean  velocities  and  h  n  is  the  mean  water  elevation.  Using 

these  substitutions  and  neglecting  any  0(d2)  terms  results  in  the  linearized  form  of  the 

shallow  water  equations  (see  Appendix  B  for  details): 

d tu*  +  Udxu  +  Vdyu*  -  f(V  +  v*)  =  -g  dxh 
dtv*  +  Udxv*  +  VdyV*  +  f(U  +  u*)  =  -g  dyh  (11.21) 

dth  +  U dxh  +  Vdyh  +  hB  (dxu*  +  dyv*)  =  0. 


F.  KLEIN-GORDON  EQUIVALENT  TO  THE  SHALLOW  WATER  MODEL 

Using  the  linearized  form  of  the  shallow  water  equations,  we  can  find  a  Klein- 
Gordon  equation  equivalent  through  a  series  of  linear  operations  as  outlined  by  Pedlosky  [13]. 
We  begin  by  defining  the  operator  (linearized  Lagrangian  derivative) 

^7  =  dt  +  Udx  +  Vdy.  (11.22) 
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Substituting  (11.22)  into  (11.21)  we  have  the  modified  form 


Dh 

~Dt 


Du* 

Dt 

=  -g  dxh 

(11.23) 

Dv* 

Dt  +/<C,+“> 

=  —  g  dyh 

(11.24) 

+  Iib  ( dxu  +  dyV  ) 

=  0. 

(11.25) 

Taking  /  as  constant  (along  some  latitude)  and  summing  the  partial  derivative  of  (11.23) 
with  respect  to  x  and  the  partial  derivative  of  (11.24)  with  respect  to  y,  we  have 


o  iDu*  \  ff)  * 
Or  (  )  -  fdxv 

'  Dv* 

d" 1  dT  1  +  fdyU * 


99xxh 

Q^yyll 


( dxu *  +  dyv*)  +  f  idyll*  -  dxv*)  =  —g'V2h 


(11.26) 


Similarly,  we  find  the  difference  of  the  partial  derivative  of  (11.23)  with  respect  to  y  and  the 
partial  derivative  of  (11.24)  with  respect  to  x  to  find 

9y  (iw)  ~~  f9yV*  =  ~gdxyh 

-  dx  +  fd*u*  =  -a^xyh 

^  (dyu *  -  dxv *)  -  /  (dxu*  +  dyv*)  =  0  (11.27) 

We  apply  the  operator  jd  to  (11.26)  and  add  — /  times  (11.27) 

^2  (d*U*  +  dVV*)  +  (dVU*  ~  d*V*)  =  ~9^  (V2^) 

+  / 2  ( dxu *  +  dyv *)  -  f-jyt  (9yU*  -  dxv*)  =  0 

^  (dxy*  +  dyv*)  +  f  (dxu*  +  dyv *)  =  -g^-t  ( V2h )  .  (11.28) 
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From  (11.25),  we  know  that  —  =  dxu*  +  dyv*,  which  can  be  used  in  (11.28) 


which,  after  factoring  f,  becomes 

wt(^+f2h-^v2h)  =  °- 

We  can  rewrite  this  equation  as 

+  fh  -  CqV2/?.  =  S(x,  y,  t ) 

where  Cq  =  ghs  and  f(^S(x,y,t)^j  =  dtS  +  UdxS  +  V dyS  =  0.  This  source  term  we 
will  assume  to  be  zero,  giving  us  the  homogeneous  form 

^  +  fh  -  c20V2h  =  0. 

If  we  expand  the  operator  twice,  we  get  an  expanded  Klein-Gordon  form 

(dt  +  Udx  +  Vdyf  h  -  CgV2/i  +  fh  =  0  (11.29) 

of  the  shallow  water  equations  under  constant  advection  U  and  V  and  dispersion  evidenced 
by  the  f  term.  This  equation  specifies  the  perturbation  of  the  wave  height  h  above  a 
reference  level  hn- 

G.  RECOVERING  THE  FLUID  VELOCITIES 

Suppose  now  that  we  have  the  solution  for  h(x,  y,  t ).  In  order  to  recover  the  fluid 
velocities  u(x,y,t )  and  v(x,y,t),  we  consider  the  modified  form  of  our  shallow  water 
model  shown  in  (II.23)-(II.25).  We  first  apply  the  operator  f  to  (11.23)  and  multiply  (11.24) 
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by  /  to  yield: 


DV  Dv*  DV  _  D  (dh 
~dT  ~  ’Jjt  ~  *~Dt  ~  ~9  Dt 
Dv*  2  dh 

S~Dt+f  {U  +  U  )_“9/Sj 


(11.30) 


(11.31) 


Adding  (11.30)  and  (11.31)  yields: 


wt+f2)u'  +  f2u 


f  D  / dh\  dh 

^  \  Dt  \  dx )  dy 


(11.32) 


The  solution  of  this  partial  differential  equation  (no  more  difficult  to  solve  than  the  equation 
for  the  perturbation  of  the  wave  height  h )  gives  us  the  fluid  velocity  in  the  x  direction. 


A  similar  manipulation  (subtracting  /  times  (11.23)  from  jh  operating  on  (11.24)) 


yields: 


4!  +  />y  +  /v^MV/ty 


Dt  \dy 


(11.33) 


The  solution  of  this  partial  differential  equation  gives  us  the  fluid  velocity  in  the  y  direction. 
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III.  HIGH  ORDER  NON-REFLECTING  BOUNDARIES 


The  numerical  solution  of  a  wave  propagation  problem  in  a  very  large  or  unbounded 
domain  provides  a  challenging  computational  difficulty  -  namely,  solving  the  problem  on 
a  finite  computational  domain  while  maintaining  the  true  essence  of  the  solution.  One  of 
the  modem  techniques  that  has  garnered  a  significant  amount  of  attention  in  handling  this 
challenge  is  the  absorbing  or  non-reflecting  boundary  condition  (NRBC)  method.  In  using 
this  method,  the  original  infinite  domain  is  truncated  by  an  artificial  boundary  B,  resulting 
in  a  finite  computational  domain  O  and  the  residual  domain  I).  Figure  4  illustrates  the 
NRBC  set-up  using  an  infinite  wave  guide.  Here,  the  artificial  boundary  B  extends  from 
the  southern  (T $ )  to  the  northern  (T jv)  boundaries  of  the  wave-guide,  thus  creating  the  east 
(r E)  and  west  (T w)  boundaries  of  Q  at  x  =  xE,  xw  respectively.  Appropriate  boundary 
conditions  are  prescribed  on  the  northern  (Tat)  and  southern  (T s)  boundaries.  Outside  of 
the  area  enclosed  by  these  boundaries  is  the  residual  infinite  domain  D. 

y 


X 

Figure  4:  An  infinite  wave-guide  truncated  by  artificial  boundaries  T w  and  T E 

One  would  expect  that  the  introduction  of  any  boundary  B,  where  one  does  not 
physically  exist,  to  pollute  the  solution  through  the  reflection  induced  by  such  an  artificial 
boundary.  In  the  last  two  decades,  significant  efforts  have  been  extended  to  find  stable, 
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efficient,  accurate  and  practical  means  of  reducing  this  reflection  through  so-called  NRBCs 
[18]. 

Several  high-order  NRBCs  have  been  devised  to  reduce  spurious  reflections  that 
would  pollute  the  solution.  Beginning  in  the  late  1980s,  the  well-known  Engquist-Majda 
[19]  and  Bayliss-Turkel  conditions  [1]  gave  way  to  Collino’s  [2]  low  derivative,  auxiliary 
variable  formulation  for  the  2D  scalar  wave  equation.  This  sparked  a  flurry  of  activity  in  an 
effort  to  find  quality,  high-order  NRBCs  that  were  easily  implementable.  The  sheer  volume 
of  literature  on  the  topic  of  boundary  conditions  for  infinite  problems  suggests  that  there 
is  no  “perfect”  boundary  condition  available  for  general  purpose  use.  In  reality,  a  modeler 
must  make  decisions  on  how  to  balance  accuracy,  efficiency  and  ease  of  implementation  to 
yield  reasonable  solutions.  Extensive  reviews  on  the  topic  can  be  found  in  [18,  20,  21,  17] 

A.  HIGDON  SCHEME 

The  starting  point  for  the  family  of  NRBCs  discussed  in  this  dissertation  is  the 
condition  devised  by  Higdon  in  a  series  of  papers  [22,  23,  24,  25,  26,  27],  that  was  demon¬ 
strated  in  a  low-order  finite  difference  setting.  While  in  theory,  Higdon’s  NRBC  is  con¬ 
sidered  a  high-order  NRBC,  the  formulation  requires  evaluation  of  increasing  high-order 
spatial  and  temporal  derivatives  as  the  order  of  the  NRBC  is  increased.  Higdon’s  condi¬ 
tion  (and  most  NRBCs  for  that  matter)  seeks  to  annihilate  waves  impinging  normal  to  a 
boundary.  To  see  the  idea  behind  this  condition,  consider  a  one-dimensional  wave  equation 

dtth  Cq  dxxh  0 

whose  solution  was  obtained  by  d’Alembert  in  1747  [28],  as 

h(x ,  t )  =  F(x  —  c0t)  +  G(x  +  c0t). 

This  solution  implies  that  there  are  two  components  to  the  solution  -  one  wave  G  of  fixed 
shape  moving  to  the  left  at  velocity  —  cq  and  one  wave  F  of  fixed  shape  moving  to  the 


24 


right  at  velocity  c0.  Now,  suppose  that  the  right  moving  wave  approaches  a  boundary.  To 
perfectly  absorb  the  wave  impinging  on  the  boundary,  the  boundary  must  satisfy  G  =  0 
such  that  the  boundary  condition  is  h(x,t )  =  F(x  —  c^t).  Differentiating  the  boundary 
condition  with  respect  to  x  and  t  results  in: 

dxh  =  F'{x  —  c0f),  dth  =  — CqF'(x  —  c0t),  (III.  1) 


which  implies 


dth  +  c0dxh  =  0. 


(HI-2) 


This  is  the  Sommerfeld  radiation  condition  for  the  eastern  boundary.  If  we  expand  the 
discussion  to  two-dimensional  problems,  this  condition  implicitly  assumes  that  by  the  time 
the  wave  front  reaches  the  boundary,  it  is  traveling  primarily  as  a  plane  wave  at  speed  c0. 


1.  Accounting  for  Dispersion 

In  a  non-dispersive  medium,  waves  can  propagate  without  deformation.  The  chal¬ 
lenge  associated  with  dispersive  waves  such  as  the  Klein-Gordon  equivalent  under  consid¬ 
eration  here,  is  that  the  wave  speed  depends  on  the  frequency  of  the  wave.  Thus,  if  using  the 
Sommerfeld  radiation  condition  (III. 2),  only  the  waves  traveling  at  phase  speed  c0  will  be 
absorbed  -  for  all  others,  only  a  portion  of  the  wave  will  be  absorbed.  Higdon  considered 
a  composition  of  J  of  these  simple  first  order  operators  to  yield  his  boundary  condition: 


Hj: 


h  =  0 


on  T 


(HL3) 


where  dn  is  the  normal  derivative  on  the  boundary.  In  the  case  of  the  wave  guide  shown  in 
Figure  4,  this  derivative  is  dx  and  —dx  corresponding  to  the  eastern  and  western  boundaries 
respectively.  The  boundary  condition  contains  parameters  Cj  that  can  be  interpreted  in 
terms  of  the  phase  velocities  of  waves  absorbed  exactly  at  the  boundary.  Except  in  contrived 
examples,  there  are  infinitely  many  waves  composing  the  solution,  and  in  a  dispersive 
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medium,  a  corresponding  infinitely  many  phase  velocities.  A  choice  of  the  order  J  of  the 
boundary  condition  seeks  to  annihilate  the  “most  significant”  J  waves. 

2.  Reflection  Caused  by  the  Boundary 

For  purposes  of  illustration,  consider  a  simplified  version  of  (11.29)  where  U,  V  are 
both  zero 

dtth  -  c20V2h  +  f2h  =  0.  (III. 4) 

Further,  suppose  that  the  domain  is  structured  such  that  the  NRBC  is  imposed  on  only  on 
the  eastern  boundary  of  a  rectangular  grid  as  shown  in  Figure  5.  On  the  south  and  north 


rN 


Figure  5:  A  semi-infinite  channel  truncated  by  artificial  boundary  T e 
boundaries  (r^  and  T N),  we  have  no  normal  flow,  i.e.,: 

dyh  =  0  on  Tat  and  (III. 5) 

and  we  impose  h  =  f(y,  t)  on  Tw-  As  x  — >  oo  the  solution  is  known  to  be  bounded  and 
not  to  include  any  incoming  waves.  A  solution  to  (III.4)  has  the  form 

h(x,  y,  t )  =  Y(y)  cos  ( kx  —  ut  +  4>)  (III. 6) 

such  that  Y(y)  satisfies  (111.4).  One  such  example  Y(y)  =  cos  ^  for  n  —  0,1,2,...  that 
satisfies  these  boundary  conditions  is  given  by  Givoli  and  Neta  [29].  Given  this  choice  for 
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Y(y),  the  dispersion  relation  is 


w2  =  cl  (V  +  +  f2,  n  =  0, 1,  2, ... .  (IIL7) 

In  this  solution,  k  is  the  horizontal  wavenumber,  n  is  a  parameter  for  controlling  the  shape 
of  Y(y),  lo  is  the  wave  frequency  and  0  is  the  phase  shift.  The  horizontal  phase  velocity  [28] 
is  therefore  Ck  =  f  for  a  particular  wave  number.  Suppose  that  one  of  the  G'/s  in  (III. 3) 
equals  Ck. 


dth  =  c oY(y)  sin  (kx  —  t ot  +  (ft) 

dxh  =  — kY(y )  sin  (kx  —  cut  +  (ft) 
to 

so  that  0  =  dth  +  —  dxh 

k 

thus,  satisfying  the  Higdon  boundary  condition  (III. 3)  exactly  for  that  particular  mode. 
If,  however,  none  of  the  Cj’s  were  identically  Ck,  then  a  portion  of  the  mode  would  be 
reflected  and  the  boundary  condition  would  not  be  exactly  satisfied.  To  make  the  boundary 
condition  true,  it  would  have  to  be  adjusted  to  incorporate  the  reflected  modes.  Higdon  [27] 
sketches  and  Dea  [17]  details  via  a  simple  induction  argument  that  proves  the  reflection 
coefficient  Rj  is 

*=ri 

3= 1 

We  notice  here  that  Rj  is  a  product  of  factors  that  are  each  less  than  one.  Therefore, 
simply  increasing  the  order  of  the  NRBC  (J)  reduces  the  amplitude  of  the  reflected  wave 
irrespective  of  the  choice  of  Cj.  van  Joolen  et  al.  [30,  page  1045]  notes,  “Of  course,  a 
good  choice  for  the  Cj  would  lead  to  better  accuracy  with  a  lower  order  J,  but  even  if  the 
‘wrong’  C/s  are  taken  . .  .one  is  still  guaranteed  to  reduce  the  spurious  reflection  as  the 
order  J  increases.”  We  can  use  the  dispersion  relation  together  with  the  definition  of  the 


Cj  -  Ck 


Ci  +  ck 


(HI- 8) 
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horizontal  phase  velocity  Ck  to  find  that 


Ck 


+ W) + /2 

k 


r2  n2  7T2  ,  f2 

L0  b2  '  J 

k2 


which  shows  that  Ck  >  c0  and  therefore  guides  our  selection  of  C3  to  always  be  at  least  c0. 

Based  on  the  preceding  discussion,  we  outline  some  of  the  inviting  characteristics 
of  the  Higdon  NRBC. 

•  The  boundary  condition  is  tractable,  extending  basic  principles  to  arrive  at  the  final 
condition. 

•  They  are  exact  for  all  waves  that  propagate  with  speed  Cr 

•  Reflection  is  guaranteed  to  decrease  by  simply  increasing  the  order  of  the  NRBC. 

•  They  have  been  applied  to  a  wide  variety  of  wave-type  problems  including  those  in 
dispersive  media  [17,  22,  23,  24,  25,  26,  27,  29,  31,  32]. 

The  boundary  condition,  however,  suffers  from  an  implementation  point  of  view.  Due  to 
the  increasing  order  in  the  spatial  and  temporal  derivatives  required,  even  Higdon  in  his 
original  papers  considered  practical  implementation  to  be  no  more  than  J  =  3. 


B.  HIGDON  ADJUSTMENTS  FOR  ADVECTION 


In  the  discussion  of  Higdon’s  boundary  condition  above,  we  considered  the  zero 
advection  case.  Eventually  we  wish  to  implement  the  Klein-Gordon  equivalent  which  in¬ 
cludes  constant  advection.  Here,  we  discuss  modifications  to  Higdon’s  scheme  to  accom¬ 
modate  non-zero,  constant  advection.  This  discussion  will  use  physical  arguments  that  will 
be  reinforced  by  numerical  considerations  later  in  this  dissertation.  To  begin,  we  consider  a 
wave  augmented  by  advection  impinging  on  the  NRBC  of  the  semi-infinite  channel  shown 
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in  Figure  5.  Implicit  in  this  derivation  is  that  by  the  time  the  wave  pulse  arrives  at  the 
NRBC,  it  is  traveling  primarily  as  a  plane  wave  (in  the  x— direction).  The  wave  moves 
according  to  the  equation  under  consideration, 

(d t  +  Udx  +  Vdyf  h  -  c20V2h  +  f2h  =  0. 

If,  however,  this  wave  is  moving  primarily  in  the  a:— direction,  then  any  effects  of  y  are 
negligible.  Further,  as  in  the  derivation  of  the  Sommerfeld  radiation  condition  (III. 2),  we 
first  consider  a  non-dispersive  environment  such  that  f2  =  0  in  this  derivation  to  suggest 
an  appropriate  boundary  condition.  This  simplifies  our  discussion  to  a  wave  that  moves 
according  to 

(dt  +  U dxf  h  -  c20dxxh  =  0.  (HI. 9) 

As  outlined  in  Appendix  C,  the  solution  takes  the  form 

h(x,t )  =  F^x  —  (c0  +  U)tj  +  G^x  +  (c0  —  U)tj  (III.  10) 

with  the  interpretation  that  the  general  solution  is  the  sum  of  F,  a  wave  of  fixed  shape 
moving  to  the  right  with  velocity  c0  +  U  and  G,  a  wave  of  fixed  shape  moving  to  the  left 
with  velocity  c0  —  U. 

As  described  in  III.A,  if  we  consider  the  wave  moving  to  the  right  approaching  the 
eastern  boundary,  the  boundary  must  satisfy  G  =  0  such  that  the  boundary  condition  is 
h(x,  t )  —  f[x  —  (c0  +  U)tj .  Differentiating  the  boundary  condition  with  respect  to  x  and 
t  results  in: 

3xh  =  F'(x-  (c0  +  U)t)  dth  =  -(co  +  U)F'(x  -  (c0  +  U)tj, 

which  implies 

dth  +  (cq  +  U)dxh  =  0. 
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Following  the  discussion  for  the  reflection  coefficient,  this  suggests  the  “educated”  choice 
of  augmenting  the  Higdon  parameter  Cj  with  the  added  advection.  This  choice  will  be 
revisited  later  in  this  dissertation. 

C.  GIVOLI-NETA  AUXILIARY  VARIABLE  LORMULATION 


In  [29],  Givoli  and  Neta  directly  extended  the  Higdon  scheme  to  high-order  finite 
difference  discretizations  via  an  algorithm  where  the  order  of  the  NRBC  was  simply  an 
input  parameter.  They  later  extended  this  formulation  [33]  to  one  that  does  not  involve 
any  high  derivatives  (hereafter  referred  to  as  the  G-N  formulation).  The  elimination  of  all 
high-order  derivatives  is  enabled  through  the  introduction  of  special  auxiliary  variables  on 
B.  This  construction  demonstrated  in  [33]  and  [34]  for  finite  differences  was  further  ex¬ 
tended  in  [35]  for  finite  element  schemes  to  solve  the  dispersive  wave  equation.  Hagstrom 
and  Warburton  [36]  also  used  the  Higdon  and  auxiliary  variable  framework  to  develop  a 
symmetric  boundary  formulation  in  a  full-space  configuration  where  special  comer  com¬ 
patibility  conditions  were  developed  for  the  non-dispersive  wave  equation.  Extensions  and 
comparisons  between  the  two  methods  were  published  by  Givoli  and  Hagstrom  et  al.  in 
[37]  and  [38], 

We  present  a  brief  summary  of  the  G-N  auxiliary  variable  process  as  described  in 
[35].  For  the  semi-infinite  channel  shown  in  Figure  5,  this  auxiliary  formulation  begins 
with  the  Higdon  boundary  condition  given  by: 


Hj : 


h 


0  on  T#. 


(HI  11) 


Auxiliary  functions  d)\,  . . . ,  4>j- i,  which  are  defined  on  Tg  as  well  as  in  the  exterior  do¬ 
main  D  are  now  introduced.  Eventually,  we  shall  use  these  functions  only  on  T e,  but  the 
derivation  requires  that  they  be  defined  in  D  as  well,  or  at  least  in  a  non-vanishing  region 
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adjacent  to  T e-  The  functions  0j  are  defined  via  the  relations 


(px  +  h  —  4>\  , 

dx  +  =  (t)2  i 


(III.  12) 
(III.  13) 


(dx  +  -^-dt')<i>j-i  =  0.  (III.  14) 

By  definition,  these  relations  hold  in  D,  and  also  on  T^.  It  is  easy  to  see  that  (III. 12  - 
III.  14),  when  imposed  as  boundary  conditions  on  T e,  are  equivalent  to  the  single  boundary 
condition  (III.  11).  If  we  also  define 


4>o  =  h  <pj  =  0  , 


(III.  15) 


then  we  can  write  (III. 12  -  III. 14)  concisely  as 

(<9z  +  ^r<9, (pj-i  =  4>j  j  =  (III.  16) 

This  set  of  conditions  involves  only  first-order  derivatives.  However,  due  to  the  appearance 

of  the  ^-derivative  in  (III.  16),  one  cannot  discretize  the  <i>3  on  the  boundary  T E  alone. 

Therefore  we  shall  manipulate  (III.  16)  in  order  to  get  rid  of  the  ^-derivative. 

The  function  h  satisfies  the  dispersive,  advective  wave  equation  (11.29)  in  I).  Since 

the  function  cf>i  is  obtained  by  applying  the  linear  (constant  coefficient)  operator  {px  + 

to  h,  it  is  can  be  shown  that  0\  should  also  satisfy  the  same  equation  in  I)4.  Further,  since 

cj)j  is  obtained  by  applying  the  same  linear  operator  j  —  1  times  to  0i,  the  functions  03 

4Here  we  must  use  the  assumption  that  c(J  and  /  are  constants.  By  applying  the  differential  operator  to 
(11.29),  computing  each  of  the  0,  derivatives  present  in  (III.  17)  using  the  differential  operator  and  simplifying, 
a  simple  induction  argument  shows  that  the  <f)j’s  must  satisfy  (III.  17) 
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should  satisfy  an  equation  like  (11.29),  namely, 


(dtt  +  (U~  —  Cq)  dxx  +  —  Cq)  dyy-\- 

2  Udxt  +  2V  dyt  +  2f/Uc^  +  /2)  y  =  0  (III.  17) 

Using  the  following  identities: 

dxx4>j  —  (dx  —  T7  )  ( &X  +  ^7  J  4>j  +  ^72  4 

\  1  /  V  W+1  /  L'j+1 

()xl0j  ((/  (<9x0j) 

^xy0j  (  A/  (  ((.r  0j  ) 


and  combining  with  (III. 16)  allows  us  to  write  (III.  17)  as: 


/2f7  U2-c2\ 

U  C?  J 


4'-i  + 


2UU 


-  21")  -  (l"2  -  C2) 

-  2t/)  -  2UV tf>'j  -  Z2^-1  = 


(f/2  -  c2)  0j+1 


for  j  =  (III.  18) 


The  details  of  this  transformation  are  shown  in  Appendix  D.  In  (III.  18)  and  elsewhere,  a 
prime  indicates  differentiation  with  respect  to  y  along  T e,  be.,  the  tangential  derivative 
along  T  f  .  As  desired,  the  new  boundary  condition  (III.  18)  does  not  involve  ^-derivatives. 
In  addition,  there  are  no  high-?/  or  t  derivatives  beyond  second  order.  It  should  be  noted  that 
van  Joolen  et  al.  [31]  developed  an  equivalent  formulation  in  using  Lagrangian  derivatives. 
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Rewriting  (III.  12),  (III.  18)  and  (III.  15),  the  new  formulation  of  the  J^-order  NRBC 
onTg  can  be  summarized  as  follows: 


Poh  +  d xh  =  , 

aj4>j- 1  +  ~~  \fij- 1  +  Pj&j  —  7 ■  f2(Pj- 1  = 


=  h 


<Pj  =  0 


(III.  19) 
(III.  20) 
(III.21) 


where 


A) 

A 


5  q 


1  21/  l/2  -  <2 

—  n  ■  = - 1 - - 

r,  ’  ■'  Cj  c]  ’ 

(c/2  -  c2)  (7  +  7- j  -  2t/,  7  =  2W, 


K<  =  77  _  2V,  A„  =  l/2  -  4 


A7.  =  U2  -cl 


What  remains  is  to  link  the  boundary  condition  to  the  interior  formulation.  As  will  be 
shown  in  the  next  chapter,  the  means  to  link  the  two  formulations  will  naturally  follow 
from  the  spectral  element  framework. 
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IV.  DISCRETIZATION  VIA  SPECTRAL  ELEMENTS 


Once  a  suitable  NRBC  is  devised,  the  problem  must  be  solved  numerically  in  O  by 
the  finite  difference,  or,  as  in  the  case  of  this  analysis,  the  spectral  element  (SE)  method. 
The  SE  method,  originally  introduced  by  Patera,  . .  combines  the  generality  of  the  finite 
element  method  with  the  accuracy  of  spectral  techniques  ...”  [39,  p.  468].  The  key  to  the 
spectral  element  (SE)  method  is  the  careful  selection  of  the  integration  and  interpolation 
points  in  order  to  yield  accurate  but  efficient  solutions. 

As  indicated  in  Chapter  III,  the  Givoli-Neta  auxiliary  variable  formulation  has  been 
previously  demonstrated  in  both  finite  difference  and  finite  element  schemes  to  arbitrarily 
high  NRBC  order,  however,  accuracy  gains  realized  by  increasing  the  NRBC  order  slowed 
significantly  after  certain  orders.  The  SE  formulation  used  in  this  dissertation  seeks  to  rem¬ 
edy  this  limitation  by  using  a  high-order  treatment  of  space  (SE)  to  show  the  benefits  of 
using  a  high-order  boundary  (G-N)  scheme.  It  should  be  noted  that  the  only  other  spec¬ 
tral  element,  high-order  boundary  approach  (to  the  author’s  knowledge)  is  [40]  that  is  in 
press.  That  paper  shows  how  a  SE  approach  combined  with  high-order  boundary  treatment 
significantly  improves  the  accuracy  for  the  non-dispersive  wave  equation  on  a  semi-infinite 
channel  using  the  Hagstrom-Warburton  formulation  [36].  The  key  difference  in  our  work  is 
that  we  use  high  order  space,  boundary  and  time  integration  (discussed  later  in  this  disser¬ 
tation)  in  both  a  dispersive  and  non-dispersive  wave  equation  setting.  We  show  the  details 
of  early  results  (absent  of  advection)  in  [41]. 

What  follows  in  this  chapter  is  a  brief  overview  of  the  SE  method  as  presented  by 
Giraldo  [42].  Additional  treatment  of  this  method  can  be  found  in  Giraldo-Restelli  [43]  or 
Pozrikidis  [44].  For  this  problem,  we  will  discuss  two  formulations  -  one  for  the  interior 
and  one  for  implementing  the  boundary  conditions. 
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A. 


INTERIOR  WEAK  INTEGRAL  FORMULATION 


Finite  difference  schemes  seek  to  discretize  differential  operators  based  on  Taylor 
series  representations  of  the  function  near  a  point.  Spectral  element  schemes  seek  to  solve 
the  equations  in  a  weak  integral  form.  It  can  be  shown  that  if  the  governing  equations 
satisfy  the  weak  integral  form,  then  they  also  satisfy  the  differential  form  provided  the  test 
functions  used  are  sufficiently  differentiable. 

The  weak  form  of  the  problem  in  Q  is  now  constructed  for  the  semi-infinite  channel 
described  in  III.A.2.  Following  a  standard  Galerkin  approach,  the  solution  is  sought  in  the 
space  of  test  functions, 

V  =  {h\h  E  H\ Q)  and  h  =  f(y,t )  on  Tw}.  (IV.l) 


Here,  is  the  Sobolev  space  of  functions  whose  first  derivatives  are  also 

square-integrable  in  O.  Now,  Equation  (11.29)  is  multiplied  by  the  globally  defined  basis 
functions  T y)  G  V  and  integrated  over  O  so  the  weak  form  is: 


f  ■■  f  d2h 

/  *ihdn  +  Xx  /  T, 

In  Jn  c,x 


(IQ  +  Xy  /  ^  i 


d2h 
dy 2 


dQ 


+2 U  I  ^~dtt  +  2V  [  ^~dn 
Jn  (xx  Jn  &  V 

+UV  I  T,  (t)'h  ■  n'h  )  dtt  +  f2  [  T,//  r/o  =  0. 

Jn  \oxoy  ay  ox )  Jq 


(IV.2) 


Here,  h  and  h  are  the  first  second  derivatives  of  h  with  respect  to  time.  Further,  ||  and 
are  the  mixed  second  derivatives  of  h  with  respect  to  t ,  x  and  y. 

In  order  to  maintain  a  symmetric  form  of  the  problem,  the  mixed  derivative  has 
been  split  appropriately.  To  ensure  the  solution  h  is  in  the  vector  space  II 1  requires  some 
special  handling  of  the  second  order  derivatives,  which  is  facilitated  by  the  use  of  Green’s 
theorem,  i.e., 


T  Q2k 

^.t7—  dU  = 

ox- 


T  dh  iti 
Xl't—nx  dF  - 
ox 


d^i  dh 
dx  dx 


dil 
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where  n  is  the  outward  normal  on  Y  and  nx  is  the  ^-component  of  that  outward  normal. 


Extending  this  idea  for  each  second  order  derivative  in  (IV.2)  gives  the  weak  form: 


f  t  v  ,  /‘  dh  ^  r  d^i  dh  ^ 

/  ^ih  dU~  Xx  /  -7—  —  dU  —  Xy  /  -77—77-  dQ 
Jn  Jn  dx  dx  Jn  dy  dy 

-uv  l  +  dn 

Jn  \  dy  ox  ox  dy  J 

+2 u  I  v~dn  +  2V  I  v—dn  +  f  I  ^thdn  (iv.3) 

Jn  dx  Jn  ay  Jn 

f  dh  f  dh  f  ( dh  dh  \ 

+AX  /  —nx  dr  +  A,,  /  ^iT^-riy  dT  +  UV  /  'Et  —ny  +  —  nx  dr  =  0. 

Jr  ox  UJ  r  dy  Jv  \dx  dy  J 


Recall  that  for  the  semi-infinite  channel  described  in  Figure  5,  the  corresponding 
boundary  conditions  for  T N  andTs  are  dyh  =  0,  displacement  is  specified  on  and 
r s  is  a  non-reflecting  boundary.  Using  these  boundary  conditions  on  the  northern  and 
southern  borders  along  with  the  normal  vectors  specified  by  the  structured,  rectangular 
grid  shown  in  Figure  6,  the  problem  may  be  simplified  to  account  for  contributions  on 
individual  boundaries.  Using  this  information  along  with  (III. 19)  to  eliminate  the  normal 
derivative  term  on  Te,  we  get: 


Figure  6:  Normal  Derivatives  to  Boundaries 
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in 


T  V  dh  . 

\Efi/r  cZfi  —  At  /  — — —  dCl  —  \v 
ox  ox 


chi',  dh 


in  u-h  u-h  Jn  dy  dy 

_uv  f  fd%dh  +  d%dh 

Jn  V  °y  dx  dx 


d,Q 

dQ 


+2 U  I  'S>~dtt  +  2V  I  '\’j>h  dil  -  f2  I  ^ihdn 

Jn  dx  Jn  dy  Jn 


+Aa 


The) i  dT  i 


dr  E 


(IV.4) 


P  0  7 

+uv  yt—drN-uv.  . 

Jrv  ^  4rs 


Vi^dTs  +  UV  [  *~drE  =  0. 
JrE  dy 


Note,  that  two  of  the  boundary  integrals  (other  than  those  on  T  E)  survive  after  simplification 
-  namely  on  T ^  and  Y s ■  This  occurs  since  the  boundary  condition  (III.5)  only  specifies  no 
flux  (i.e.,  reflecting  boundary  conditions)  on  these  boundaries.  This  is  a  common  boundary 
condition  for  inviscid  flow  problems.  5 


B.  BOUNDARY  WEAK  INTEGRAL  FORMULATION 


Since  the  term  <pi  appears  in  the  interior  formulation  (IV.4),  this  is  not  a  complete 
weak  form  of  the  problem  in  Q.  We  turn  our  attention  to  computing  a  separate  weak  form 
for  (III. 20)  to  complete  the  problem  statement.  This  solution  is  sought  in  the  space  of  test 
functions, 

VrB  =  {4>Mi  e  H\rE)}  (IV.5) 

where  f/1(rs)  is  the  Sobelov  space  of  functions  whose  first  derivatives  are  also  square 
integrable  onTB. 

As  in  the  interior  formulation,  we  multiply  (III. 20)  by  the  globally  defined,  one¬ 
dimensional  basis  functions  (on  the  boundary)  Q  e  VrE  and  integrate  it  over  T E.  After 

5It  should  be  noted  that  in  this  analysis,  the  two-dimensional  basis  functions,  'I-',,  are  constructed  using 
tensor  products  of  one-dimensional  basis  functions,  Q.  This  means  that  on  the  boundaries,  these  T,  basis 
functions  collapse  to  Q  on  the  boundary.  Therefore,  in  practice  the  boundary  integrals  are  constructed  using 
(i  basis  functions,  requiring  far  less  storage. 
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integration  by  parts  and  simplifying  this  yields: 


a,- 


Ci4>j-1  d^E  +  Kj  /  007  —  1  d^E  +  \  /  C'iftj- 1  dVE  +  f3j  /  Q(j)j  dT E 


-7/  (ifij  dr E  -  f2  drE  =  \x  / 


(iv.6) 


for  j  =  1, . . . ,  J  —  1.  As  in  the  interior  formulation  (IV.2),  0  and  0  are  the  first  and  second 
derivatives  of  0  with  respect  to  time  and  q Y  is  the  tangential  time  derivative  of  q b.  Recall 
from  the  auxiliary  variable  formulation  (III. 21)  that  0O  =  h  and  (j>j  =  0. 

The  formal  problem  statement  is  then:  Find  h  G  V  and  G  Vr  E  where  j  = 
1, . . . ,  J  —  1,  such  that  Equations  (IV. 4)  and  (IV.6)  are  satisfied  V  47  G  V  and  Q  G  VrB. 

C.  GRID  GENERATION  AND  CHOICE  OF  BASIS  FUNCTIONS 


One  of  the  key  advantages  of  the  spectral  (finite)  element  formulation  over  differ¬ 
encing  schemes  is  the  ability  to  represent  complex  geometries  by  breaking  up  the  domain 
into  simple  geometric  shapes.  Triangles  and  quadrilaterals  are  primarily  used  to  represent 
these  geometries  in  two-dimensions,  while  tetrahedra  and  hexahedra  are  primarily  used  to 
represent  geometries  in  three-dimensions.  The  sheer  volume  of  grid  generation  software 
available  for  generating  triangles  (tetrahedra)  and  the  dearth  of  grid  generation  software  for 
quadrilaterals  (hexahedra)  suggests  that  the  former  are  the  more  common  means  of  repre¬ 
senting  the  geometry  of  a  problem.  However,  there  are  some  very  nice  properties  of  the 
latter  that  led  to  the  choice  of  structured  (and  later  unstructured)  quadrilaterals  to  discretize 
the  geometry  in  this  analysis. 

While  triangular  grids  are  typically  easier  to  generate,  quadrilateral  elements  can 
be  formed  directly,  or  more  commonly,  by  combining  two  or  four  basic  triangle  elements 
as  shown  in  Figure  76.  This  idea  can  then  be  used  to  find  an  appropriate  mesh  for  simple 
or  complex  geometries.  For  this  analysis,  unstructured  quadrilaterals  were  generated  using 
6Adapted  from  [45],  Figure  5.2,  p.  141. 
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Figure  7:  Arbitrary  quadrilateral  element  formed  by  combining  triangles. 


Automesh-2D  [46],  an  automatic  mesh  generator  that  provides  high-quality  results  for  the 
domains  considered  in  this  dissertation.  A  sample  of  the  mesh  that  can  be  generated  is 
shown  in  Figure  87.  Noteworthy  in  this  mesh  is  that  the  complex  geometry  of  each  letter 
has  been  broken  up  into  approximately  200  non-overlapping  quadrilateral  elements  Qe 
such  that 

Nett  200 

o  =  (J  oe. 

e=l 

Further,  each  Qe  is  constructed  such  that  all  internal  angles  are  all  between  30°  and  150°  - 
a  key  criterion  in  providing  a  stable  numerical  solution. 


Figure  8:  Sample  mesh  generated  using  Automesh-2D. 

To  see  why  this  criterion  is  important,  consider  a  mapping  from  the  physical  co¬ 
ordinate  space  to  a  canonical  element  space  such  that  x  =  x(£,  //)  and  y  =  y(^,  rj)  where 
7Geometry  for  the  letters  graciously  provided  by  [46]. 


(£,  y)  e  Qr  =  [—1,  l]2.  This  nonsingular  mapping  is  made  to  facilitate  efficient  and  accu¬ 
rate  computation  of  operations  such  as  differentiation  and  integration  [47].  One  such  ele¬ 
ment  that  illustrates  this  mapping  can  be  seen  in  Figure  9.  All  derivatives  are  then  mapped 


y 


Figure  9:  Mapping  from  physical  space  to  computational  space 


to  the  canonical  element  space  by  virtue  of  the  chain  rule  to  yield  the  total  derivative: 

(  dx\  /  if  I5  \  /  df  \ 

UHtiJU)' 

From  the  total  derivative,  we  can  find  key  components  required  for  the  evaluation 
of  the  integrals  induced  by  this  canonical  element  mapping  [42].  Such  terms  include  the 
transformation  Jacobian  and  its  associated  determinant: 


(dx  dx  \ 

di  dr, 

dy  dy  J  ’ 

dt  dr,  J 


dx  dy  dy  dx 
dC,  dy  dC,  dy  ’ 


as  well  as  mappings  of  metric  terms: 


(IV.7) 


d £  1  dy  d£  1  dx  dy 

dx  |  Je |  dy1  dy  \Je\  dy  ’  dx 


1  dy  dy  1  dx 
Je\d£1  dy  |  Je  |  <9£  ’ 


(IV.8) 


If  any  of  the  nodes  comprising  the  elements  approach  interior  angles  of  180°  (degenerating 
the  quadrilateral  into  a  triangle),  it  can  be  shown  that  the  transformation  Jacobian  at  that 
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point  tends  to  0  causing  each  of  the  metric  terms  to  tend  toward  infinity  thus  corrupting  the 
solution.  Derivation  of  the  metric  term  mappings  and  consequences  of  degeneration  of  the 
quadrilateral  are  detailed  in  Appendix  E. 

Crucial  to  the  construction  of  SE  approximations  is  the  representation  of  the  solu¬ 
tion  variables  in  terms  of  smooth  basis  functions  While  in  theory  these  basis  functions 
are  defined  throughout  the  entire  domain  they  have  compact  support,  and  are  in  practice 
constructed  locally.  This  implies  that  we  can  solve  the  global  problem  by  simply  summing 
up  the  contributions  from  the  smaller  elemental  problems  in  a  process  known  as  direct 
stiffness  summation.  Ideally,  these  basis  functions  will  not  only  support  the  grid  chosen, 
but  will  also  have  properties  that  facilitate  the  numerical  integration  of  the  weak  formula¬ 
tion  [48]. 

For  this  analysis,  the  local  element-wise  solution  h  is  represented  by  Nth  order 
Lagrange  polynomials  in  (£,  rj)  such  that 

mn 

v)  =  ^  h  (&> 

k= 1 

where  (£fc,  rjk)  are  the  MN  =  (N  +  l)2  Legendre-Gauss-Lobatto  (LGL)  interpolation  points 
and  iff  are  the  associated  multivariate  Lagrange  polynomial  basis  functions.  The  square 
structure  of  the  canonical  element  simplifies  matters  in  that  we  can  express  iff  by  a  tensor- 
product  of  the  one-dimensional  Lagrange  polynomial  basis  functions  as 

(f ,  v)  =  ui  (0  Vj  {rj) 

where  iff  =  1,  2, . , . ,  TV  +  1  and  k  =  1,  2, ... ,  MN.  Further,  //,  and  i/j  are  the  one¬ 
dimensional  basis  functions  generated  using  the  LGL  points  in  £  and  //.  To  get  from  the 
one-dimensional  local  indices  {iff)  to  the  two-dimensional  local  index  k  requires  the  map¬ 
ping  k  =  i  +  (j  -1)(N  +  1). 
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The  choice  of  the  LGL  points  for  interpolant  generation  have  the  especially  attrac¬ 
tive  property  of  allowing  efficient  quadrature.  In  particular,  any  polynomial  of  up  to  degree 
2 N  —  1  can  be  integrated  exactly  (to  machine  precision)  by  evaluating  it  at  only  N  +  1 
points.  Additionally,  if  the  interpolation  points  are  also  used  as  the  sampling  points  in  the 
Gaussian  quadrature  rule,  it  has  the  added  benefit  of  yielding  a  diagonal  mass  matrix  as  the 
LGL-based  basis  functions  'F  satisfy  a  discrete  orthogonality  via  the  cardinality  property 
of  the  Lagrange  basis  functions.  The  effect  of  these  choices  is  to  allow  the  appropriately 
weighted  summation  of  the  basis  function  expansion  of  h  to  evaluate  all  integrals.  There¬ 
fore,  to  approximate  a  local  elemental  integral,  we  find 


h(x,  y)dQe 


h{t,-n)\je\dnc 


N+ 1 

~  u  (6) u  fai)  h  (&>  Vj)  \je  (&,  Vj) 
i,j= 1 


where  u  (£*)  and  c o  ( rjj )  are  the  quadrature  weights  associated  with  one-dimensional  LGL 
quadrature  and  is  the  canonical  element  region  of  integration.  This  process  is  described 
in  greater  detail  in  [42] . 


D.  GALERKIN  EXPANSION 


We  now  turn  our  attention  to  the  spatial  discretization.  First,  we  construct  the  Nth- 
order  approximation  of  the  variables  h  and  (j>j  using  the  same  basis  functions  used  in  the 
weak  form  IV.4  and  IV.6  as  follows: 

NP  Nb 

hN  =  ^2  ^kh\  4>jN  =  '22  J  =  1,2,...,  J-  1. 

k=l  k= 1 

Here,  Np  and  A),  are  the  number  of  points  that  O  and  T e  are  discretized  into,  repsectively. 
For  a  structured  quadrilateral  grid  with  N*  and  Nj1  elements  in  the  x  and  y  directions, 
Np  =  (N£N  +  1)  (N^N  +  1)  and  Nb  =  N^N  +  1.  Next,  we  substitute  this  basis  function 
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expansion  directly  into  the  weak  formulations,  resulting  in  the  following  interior 


Mh  +  (2 U Dx  +  2V  Dy)  h  +  (—A XLXX  —  A yLyy  +  f  2M  —  UV  ( Lyx  +  Lxy)^j  h  (IV.9) 
+XxAE(t)l  -  (30\xBEh  +  ( UVCN  -  UVCS  +  UVCE )  h  =  0 


and  boundary 


ajM%_ i  +  CjD%-!  +  (A yLb  -  f2Mb )  ^  (IV.  10) 

+f3jMb<j)j  -  -  K.Mb(f)j+i  =0,  j  =  1, . . . ,  J  -  1 


forms  of  the  problem.  Here,  M,  Dx,  Dy,  Lxx,  Lyy,  Lxy  and  Lyx  are  interior  formulation 
matrices  of  size  Np  x  Np.  Further,  BE,  CE,  Cs  and  CE  are  interior  formulation  matrices 
integrated  along  the  boundaries  of  size  Np  x  Np  while  AE  is  of  size  Np  x  Nb.  Finally, 
Mb,  Db  and  Lb  are  auxiliary  formulation  matrices  of  size  Nb  x  Nb.  These  global  matrices 
are  obtained  from  analogous  element  arrays  via  assembly,  given  by: 


Ne 


Ne 


M  =  f\  Me,  BE=/\  B%, 

e=l  e=l 

Ne  Ne 

LXX  —  f\  Lxx,  Lyy  =  Lyy, 
e=l  e=l 

Nb  Nb 

Mb  =  /\  mb,  Db  =  [\  db, 

6=1  6=1 

Ne  Ne 

cs  =  /\  C%,  CE  =  !\  C%, 


e=l 


e=l 


Ne 


Ne 


d„  =  A  D»  =  A  Dl 

e=l  e=l 

Ne  Ne 

Lxy  ~  /\  Lxy,  Lyx  =  Lyx, 

e=l  e=l 

Nb  Ne 

L"  =  A  Cn  =  A  C». 
6=1 

Ne 

Ae  —  f  \  Ae 


e=l 


e=l 
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Ne 

where  j\  is  the  direct  stiffness  summation  operator  required  by  all  continuous  Galerkin 

e=l 

methods.  The  expressions  for  the  element  arrays  are: 


Me-  = 
13 


De  ■  = 
y,v 


Te  = 

xy,ij 


dt  = 


s~ie  _ 

Us,ij  - 


Ai>j  d£le 


f  /  dA 

f  dA&ipj 


dy  dx 


vrv's  dFe 


dQe 
dflP 


'T§ 


$ 


rfrs 

■  dx  ■ 


ne  _ 

nE,ij  ~ 


Te  — 
^xXjij 


re  = 

yx,ij 


t  = 


s~ie  _ 

UE,ij  ~ 


I  Ai>j  dTe 

f  9AdA 

/  o  _  (Xu  Le 

'q  OX  uX 


dA 

dx  dy 


dflP 


IjpE 


A 


di±dTE 

dy  ' 


De  = 

x,%3 


Te  = 
yy,ij 


mrs  = 


/~ie  _ 

0 N,ij  ~ 


Ae 

A EAr  ~ 


f  A^r1  dtt 

!ne  dx 

f  dA  d'lp-j 


'ne  dy  dy 
[  i3 r  13 s  dr  e 


dfL 


!  piv 


V'dx  • 


I  rE 


'>Pii3r  dr 


where  Vte  and  Te  denote,  the  part  of  Q  and  V  associated  with  element  e.  Also,  A  and  u, 
are  the  locally  defined  basis  functions  from  which  the  global  basis  functions  (4q  and  Q)  are 
formed.  For  quadrilateral  elements  with  polynomial  order  N  as  used  in  this  analysis,  A  W'H 
be  discretized  into  (N  + 1)2  points,  and  z/*  into  (N  + 1)  points,  thus,  i,j  =  1, 2, . . .  (A^  + 1)2 
and  r,  s  —  1,  2, . . .  N  +  1. 

Now,  let: 


A  =  M-1  (A,A xBe  -  2 UDX  -  2 VDy) 

B  =  M-1  (A XLX  +  XyLy  -  f2M  +  UV  (Lxy  +  Lyx  -  CN  +  Cs  -  CE))  (IV.  11) 
C  =  -\xM~1AbE 


and  substitute  equations  (IV.  11)  into  (IV.9)  to  get  the  matrix  form  of  the  interior  problem: 


h  =  Ah  +  Mh  +  C0i 


(IV.  12) 
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If  we  examine  the  boundary  auxiliary  variable  formulation  (III.20),  we  see  that 
the  selection  of  appropriate  Cj  values  for  the  auxiliary  variables  has  not  yet  been  fully 
addressed.  As  discussed  in  III.A.2,  any  choice  of  Cj  is  guaranteed  to  reduce  spurious 
reflection  as  the  order  of  the  NRBC  (J)  increases,  however  based  on  the  discussion  in 
III.B,  the  Higdon  parameter  should  somehow  be  augmented  with  the  advection. 

Armed  with  this,  we  choose  convenient  values  for  our  C/s  that  cause  the  second 
order  in  time  (cty)  terms  to  vanish.  In  the  case  of  the  semi-infinite  channel,  on  the  eastern 
boundary  this  value  is:  Cj  —  c0  +  U.  A  physical  argument  for  this  choice  is  that  since  the 
predominant  speed  of  the  wave,  absent  the  advection  terms,  is  cq,  the  Cj  terms  “correct”  the 
boundary  formulation  to  account  for  the  advection.  This  selection  for  the  C'/s  then  makes: 

oli  —  ot-i  —  ■  ■  ■  —  OLJ- 1  =  0, 

«1  =  K2  =  ...  =  Kj- 1  =  K, 

Pi  —  P2  —  ■  ■  ■  —  Pj-1  =  P 

Now,  let: 

D  =  f2Mb  -  \yLb  (IV.  13) 

and  substitute  (IV.  13)  into  (IV.  10)  to  get  the  following  form  of  the  problem: 

PMbj> i  =  D  hVE  -nDbur  +  yDb  fa  +  A xMbfa 
k  Db  fa  +  pMbfa  =  IDfa  +  7 Db  fa  +  XxMbfa 
k  Db  fa  +  /3Mbfa  =  Bfa  +  yDb  fa  +  \xMbfa 

K  Db  fa-2  +  pMbfa_ 1  =  D  fa_2  +  yDb  fa-. 
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where  hrE  is  the  value  of  h  on  Y E.  If  we  now  collect  the  terms  on  the  left  and  right,  we  get 
the  matrix  form  of  the  problem: 

+  h  (IV.  14) 

where: 
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V.  SOLUTION  OF  THE  DYNAMIC  SYSTEM 


The  equations  described  by  (IV.  12)  and  (IV.  14)  constitute  a  coupled  system  of 
ODEs  that  must  be  solved  for  h(x,y,t).  Our  approach  uses  standard  klh  order  explicit 
Runge-Kutta  (RK)  methods  to  integrate  the  system  in  time.  Recall  that  RK  is  used  to  solve 
first  order  ODEs,  and  as  such,  the  second  order  systems  described  must  be  converted  into 
a  larger  system  of  first  order  ODEs.  Using  the  substitution  v  —  h  ,  v  —  h  yields  the  first 
order  systems: 


/  0 

h 

0  / 

h 

0 

— 

+ 

0  / 

V 

1  A 

V 

C01 

E<i>  =  F4>  +  h. 


(V.l) 


This  system  was  solved  using  a  two-stage  approach  at  each  time  step.  First,  the  auxiliary 
system  was  solved  to  find  the  component  <pi  required  for  the  main  system.  Then  the  main 
system  was  solved  to  find  h  at  the  next  time  step.  This  process  was  continued  until  the  final 
time  t  f  was  reached. 

A.  RUNGE-KUTTA  METHODS 


RK  methods  are  convenient  in  this  analysis  in  that  the  machinery  to  decrease  the 
truncation  error  is  the  same  for  2nd ,  3rd  or  even  kth  order  schemes.  While  there  are  other 
schemes  that  require  fewer  function  evaluations  and,  in  fact,  are  more  efficient  than  the 
explicit  RK  schemes  used  in  this  dissertation,  we  desired  a  method  that  could  change  the 
order  of  the  time  integration  method  simply  as  a  parameter  while  using  the  same  basic 
formulation.  RK  methods  use  a  standard  formulation  that  is  outlined  as  follows. 
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To  obtain  an  approximation  to  the  solution  of  the  initial  value  problem 


y'  =  f  (t,y(t)) ,  y{t0)  =  yo 

a  successive  approximation  yn+1  to  yn  is  given  by 

S 

yn+ 1  =  yn  +  At  y,  biki 


i= 1 


where  ki  =  /  fn  +  qAf,  yn  +  At  ^  a^kj  )  i  =  1,  2, . . .  s. 

V  j= i 

Here,  the  coefficients  a^-,  bi  and  c,  are  given  by  the  Butcher  tableau  for  the  given  RK 
scheme.  The  classical  explicit  RK-4  scheme  is  given  by 


Cl 

On 

a  12 

a  13 

a44 

0 

0 

0 

0 

0 

C2 

®21 

C'22 

«23 

«24 

1 

2 

1 

2 

0 

0 

0 

C3 

«31 

«32 

«33 

034 

=*•  1 

0 

1 

2 

0 

0 

c4 

04i 

a42 

a43 

Q44 

1 

0 

0 

1 

0 

bi 

b2 

h 

64 

1 

6 

1 

3 

1 

3 

1 

6 

For  this  analysis,  all  coefficients  were  computed  using  Maple™  using  routines  coded  by 
Stone  [49].  These  routines  allow  for  high  precision  computation  of  the  coefficients  re¬ 
quired  for  each  high-order  tableau.  Of  course,  the  computational  complexity  increases  for 
increased  RK  order.  To  be  specific,  the  required  number  of  function  evaluations  (stages) 
and  non-zero  parameters  required  for  the  kth  order  RK  method  are  given  in  Table  1 

Table  1:  Computational  cost  for  RK  methods  used. 


RK  Order 

2 

3 

4 

5 

6 

7 

8 

9 

10 

Stages 

2 

3 

4 

6 

7 

9 

11 

16 

18 

Parameters 

3 

8 

10 

22 

33 

45 

59 

102 

113 
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B.  SIMPLE  EXAMPLE  OF  CONVERGENCE 


To  see  the  effect  of  high-order  time  integration,  consider  the  simple  initial  value 
problem 

^  =  f{t,y)  =  (2-t)y,  y{  0)  =  e“2  (V.2) 

whose  solution  is 

y(t)  =  e-0-5^2)2. 


Now,  we  integrate  (V.2)  in  time  using  our  RK  schemes  (RK-2  -  RK-10)  from  t  —  0  to  t  —  5 
and  quantify  the  errors  observed  over  time  between  the  exact  and  RK  solution  using  the 
normalized  L2  error  defined  as  follows: 


error  \  \  ^2 


/ 


N 


n=  1 


^2  \  Vn  -  y(nAt  +  t0) 


2\ 


1 

2 


^y(nAt  +  t0)'2 

n— 1 


where  yn  and  y(nAt  +  t0)  are  the  numerical  and  reference  solutions  at  a  given  time. 

Using  a  time-step  At  =  0.05,  we  see  exponential  convergence  (to  near  machine 
precision)  of  the  error  as  indicated  in  Figure  10.  It  should  be  noted  that  in  this  example, 
even  if  an  extremely  fine  time  step  is  chosen,  the  computed  error  using  lower  order  RK 
methods  cannot  achieve  machine  precision. 

Take  for  instance  the  RK-4  method  that  has  a  truncation  error  of  O  (At4).  A  step 
size  of  At  =  10~4  should  result  in  an  error  O  (10-16),  yet  experimentally,  the  error  is 
O  (10-12).  This  apparent  discrepancy  can  be  explained  by  the  round-off  error  due  to  the 
finite-precision  arithmetic  (double)  used  in  the  computation.  This  error  increases  with  the 
total  number  of  integration  steps  used.  This  implies  that  to  integrate  (V.2)  from  t  —  0 
to  t  —  5,  we  must  evaluate  f(t,y)  a  total  of  200,000  times  for  RK-4  at  a  time-step  of 
At  =  10-4.  If  we  compare  this  with  a  time-step  of  At  =  0.025  using  RK-10,  this  should 
result  in  an  error  O  (10-16).  To  integrate  (V.2)  using  RK-10  for  the  same  period  of  time, 
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Figure  10:  L 2  error  in  solution  of  (V.2)  using  RK-2  -  RK-10  time  integration. 


we  must  evaluate  f(t,  y )  only  3,600  times.  Experimentally,  RK-10  achieves  this  O  (10~16) 
error  convergence,  thus  highlighting  the  importance  of  round-off  error  in  the  convergence 
of  the  time  integration  scheme. 

To  check  the  convergence  of  these  RK  methods  for  reasonable  time-step  values 
At  G  [0.001, 1],  consider  the  rate  of  convergence  adapted  from  [50]  and  defined  as  follows: 


rate  = 


log  (error  Atn+1)  -  log  (error Atn) 
log  (Atn+i)  -  log  (Atn) 


(V.3) 


This  rate  is  a  measure  of  how  rapidly  a  given  time  integration  method  converges  as  a  func¬ 
tion  of  the  time-step  refinement.  Figure  1 1  shows  the  normalized  L2  error  as  a  function 
of  the  time-step  refinement.  For  each  RK  order,  the  rate  of  convergence  is  averaged  over 
all  At  refinement  levels.  This  average  rate  is  annotated  for  each  RK  order  and  shows,  as 
expected,  that  the  maximum  rate  of  convergence  is  the  RK  order.  From  this  figure,  we 
see  that  this  theoretical  rate  of  convergence  is  nearly  reached  in  every  case.  It  should  be 
noted  that  once  errors  reach  O  (10-15),  round  off  errors  prevent  further  improvements  us¬ 
ing  a  given  RK  method,  and  as  already  discussed,  these  errors  actually  prevent  theoretical 
improvements  of  lower  order  (RK-2  and  RK-4)  methods. 
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Figure  11:  L2  error  as  a  function  of  time-step  refinement.  Rates  as  defined  by  (V.3)  should 
correspond  to  RK  order. 


Of  course,  one  might  wonder  why  such  high  precision  when  integrating  in  time 
would  ever  be  required.  For  most  engineering  applications,  where  the  model  parameters 
themselves  are  approximated,  this  argument  is  a  strong  one.  Perhaps  a  lower  order  time 
integrator  would  result  in  errors  that  are  on  the  same  order  as  the  model  parameters  and 
any  computational  effort  incurred  in  using  a  high-order  time  integrator  would  be  wasted. 
However,  since  this  work  centers  on  the  use  of  high-order  spectral  elements  along  with 
high-order  boundary  conditions;  high-order  time-integrators  were  also  explored  to  examine 
all  of  the  limiting  factors  of  high- order  accuracy. 

We  propose  that  in  order  to  see  the  full  improvements  of  high-order  boundary  con¬ 
ditions  requires  a  balance  of  truncation  errors  (and  round-off  errors)  between  all  of  the 
components  of  the  numerical  model;  this  includes  the  boundary  conditions,  the  spatial 
discretization  method,  and  the  time-integrators.  Practically  speaking,  we  believe  that  the 
mathematical  formulation  should  be  the  strong  point  of  any  model  evaluation,  specifically, 


53 


the  order  of  the  method  should  be  chosen  to  match  (or  exceed)  other  errors  inherent  in  the 


model.  Experiments  in  this  dissertation  were  conducted  using  boundary  conditions  up  to 
order  20,  SE  polynomials  up  to  order  16,  and  time-integrators  up  to  order  10. 
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VI.  NUMERICAL  EXPERIMENTS 


Several  numerical  experiments  were  performed  to  solve  the  Klein-Gordon  equiv¬ 
alent  (11.29)  with  and  without  dispersion  and  advection.  Additionally,  we  consider  five 
primary  configurations:  semi-infinite  channel,  infinite  channel,  quarter  plane,  open  plane, 
and  open  unstructured  plane.  Each  of  these  configurations  will  be  described  in  detail  later 
in  this  chapter. 

In  order  to  simplify  the  numerical  simulation  of  the  problem  at  hand,  the  KGE 
(11.29)  is  converted  to  anon-dimensional  form  as  described  in  [51].  Using  typical  mesoscales 
in  the  ocean,  the  length  scales  were  chosen  0(100  km),  vertical  depth  scales  0(100  m)  and 
the  dispersion  parameter  /  for  Coriolis  O(10-4s-1).  Majda  [51,  page  61]  describes  typical 
advective  velocities  as  roughly  yjk  that  of  the  medium  wave  speed  c0.  Here,  we  choose 
to  challenge  the  boundary  further  by  choosing  faster  advective  velocities  at  ^  that  of  the 
medium  wave  speed.  The  derivation  details  of  the  non-dimensional  formulation  are  given 
in  Appendix  F. 

For  experiments  that  follow,  the  reference  wave  speed  is  c0  =  1,  which  allows  the 
initial  wave  to  propagate  through  the  region  of  interest  in  a  reasonable  time  for  all  exper¬ 
iments.  Given  the  scale  choices  above  yields  a  dispersion  parameter  f2  of  0.1,  however, 
to  ensure  that  dispersion  is  felt,  we  choose  f  2  =  0.5  corresponding  to  more  than  doubling 
of  the  angular  velocity  of  the  earth.  Further,  U  and  V  are  set  to  0.0  or  0.1  under  a  two- 
dimensional  Gaussian  initial  condition  to  test  the  formulation  in  a  semi-infinite  and  infinite 
channel  setting. 

The  experiments  begin  with  an  analytic  benchmark  where  a  solution  is  synthesized 
to  the  KGE  and  compared  to  results  obtained  using  the  NRBCs.  The  chapter  continues  with 
more  general  experiments,  which  do  not  have  analytic  solutions.  In  channel  experiments, 
two  different  initial  conditions  were  tested.  The  first  is  a  two-dimensional  cosine  pulse 
where  b  is  the  height  of  the  channel  (in  all  cases  we  used  b  =  4).  This  pulse  is  chosen  with 
compact  support  such  that  the  waves  are  generated  only  in  a  narrow  region  of  the  domain 
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and  is  zero  elsewhere.  Additionally,  the  initial  condition  (by  design)  satisfies  the  no-flux 
boundary  conditions  on  the  northern  and  southern  boundaries.  The  initial  condition  is  given 
by 

h(x,  y,  0)  =  e-10*2  cos  ,  h(x,  V,  °)  =  (VI.  1) 

The  other  initial  condition  used  in  the  channel  and  remaining  experiments  is  a  two-dimensional 
Gaussian  centered  at  x  =  0,  y  =  0,  given  by 

h(x,  y,  0)  =  e-10(x2+y2) ,  h(x,  y,  0)  =  0.  (VI.2) 

In  order  to  see  the  effect  of  the  NRBC,  we  compare  our  solutions  to  the  one  com¬ 
puted  on  a  larger  domain  (for  all  experiments  except  the  analytic  benchmark).  We  consider 
this  reference  solution  “exact”  when  using  equal  order  basis  functions  on  a  fine  mesh, 
defined  in  greater  detail  in  each  of  the  subsequent  experiments.  Time  integration  was  per¬ 
formed  using  various  order  Runge-Kutta  schemes  using  a  time-step  chosen  to  ensure  a 
Courant  number  of  0.25,  where  the  Courant  number  is  conservatively  defined: 

Cn  At 

Courant  number  =  — ;  (VI. 3) 

\j  (Ax)2  +  (Ay)2 

Here,  Ax  and  Ay  are  chosen  as  the  minimum  distance  between  any  two  points  in  the 
x  or  y  directions  respectively.  This  choice  is  made  since  the  interpolation  points  are  not 
uniformly  distributed  when  using  spectral  elements.  Consider  a  canonical  element  using 
8th  order  basis  functions  with  points  chosen  using  the  LGL  interpolation  points  described 
in  Chapter  IV.  As  shown  in  Figure  12,  these  points  are  distributed  in  such  a  way  to  facilitate 
interpolation  and  integration  and  are  not  uniformly  distributed. 
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Figure  12:  Canonical  element  using  8th  order  basis  functions  showing  distribution  of  grid 
points. 


In  order  to  quantify  the  errors  observed  between  the  reference  and  NRBC  solutions, 
we  use  the  normalized  lj}  error  defined  as  follows: 


error  r 2  = 
1  UL,n 


(VI.4) 


where  Np  is  number  of  points  in  Q  and  h  P  and  are  the  numerical  and  reference  solutions 
at  point  k.  What  follows  is  a  series  of  experiments  that  were  designed  to  demonstrate  the 
efficacy  of  the  G-N  NRBCs. 


A.  ANALYTIC  BENCHMARK  SOLUTION  OF  SEMI-INFINITE  HORIZON¬ 
TAL  CHANNEL 


Recall  the  semi-infinite  channel  formulation  discussed  in  Chapter  IV  that  we  now 
implement.  Specifically,  Yw  is  introduced  at  x\y  =  —2  and  no-flux  boundary  conditions 
are  specified  on  T ^  and  T s  such  that  dyh  =  0.  Finally,  the  G-N  NRBC  B  is  introduced  at 
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xe  =  2  and  advection  in  the  x  and  y  directions  are  set  to  zero.  This  specific  situation  is 
shown  in  Figure  13. 


Figure  13:  The  semi-infinite  channel  domain  under  consideration.  Domain  is  truncated  by 
an  artificial  boundary  B  at  xe  =  2. 


A  priori,  we  synthesize  a  solution  that  satisfies  the  KGE  with  zero  advection.  Kucherov 
and  Givoli  [40]  use  a  similar  benchmark  in  a  non-dispersive  environment  for  a  single  wave 
mode.  The  solution  used  here  is  a  linear  combination  of  two  waves  and  has  the  form 

2 

hbm(x,  y,t)  =  cos  cos  [kmX  _  Umt  +  0m)  (VI. 5) 

m= 1 

where  the  parameters  are  the  same  as  defined  in  (111.6).  Given  choices  for  k,  b,  c0,  n,  /  and 
0,  we  can  determine  u  to  satisfy  the  KGE.  Here,  we  choose  /c  =  {|,|},6  =  4,  Cq  =  1, 
n  =  (2, 4}  ,  f2  =  0.5  and  0  =  (0,  zr}.  These  choices  ensure  that  the  no-flux  boundary 
conditions  on  T ^  and  T s  are  upheld  and  the  gradients  on  either  side  of  the  NRBC  are 
matched  at  t  —  0.  For  the  NRBC  solution,  we  specify  the  values  for  hrw  based  on  the 
analytic  solution.  These  parameter  choices  result  in  horizontal  phase  velocities  C0  of  1.48 
and  4.22  for  each  of  the  waves. 

The  NRBC  parameters  are  selected  simply  as  C3  =  c0,  Vj  as  described  in  Chap¬ 
ter  IV. D.  Recall  that  this  choice  eliminates  the  second  order  time  integration  terms  in  the 
boundary  formulation.  In  theory,  if  the  C3  parameters  were  chosen  to  match  the  horizon- 
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tal  phase  velocities,  J  =  2  would  be  highest  order  NRBC  required.  In  general  problems, 
however,  these  phase  velocities  are  not  normally  known  a  priori.  We  therefore  keep  our  for¬ 
mulation  general,  relying  on  the  fact  that  simply  increasing  the  NRBC  order  is  guaranteed 
to  reduce  the  reflection  caused  by  the  boundary. 

1.  Results 

In  Figure  14,  we  show  the  reference  solution  on  the  top  panel  and  the  solution  on 
the  J  =  4  NRBC  truncated  domain  on  the  bottom  panel  at  t  =  3.  The  NRBC  solution  uses 
4th  order  spectral  elements  on  a  24  x  24  element  grid,  discretizing  the  domain  into  9,409 
points.  Qualitatively  speaking,  the  results  show  very  little  reflection  when  compared  to  the 
synthesized  solution. 

F 

0 

-0.5 

I::, 


Figure  14:  Semi-infinite  channel  comparing  synthesized  solution  and  the  NRBC  solution. 
4th  order  spectral  elements  using  NRBC  order  J  =  4  with  zero  advection  at  t  —  3  are 
shown. 


Synthesized  Solution 


X 


NRBC  Solution  (J=4) 


Quantitative  results  can  be  observed  in  Figure  15  showing  the  error  on  as  a  func¬ 
tion  of  spectral  and  NRBC  order  (J  =  1, . . . ,  10, 15,  20).  The  number  of  elements  is  ad¬ 
justed  for  each  spectral  order  to  maintain  an  equal  number  of  points  (9,409)  that  the  domain 
is  discretized  into.  It  is  clear  that  increasing  the  NRBC  order  yields  significant  gains  in  ac¬ 
curacy,  but  by  J  =  5,  the  spatial  discretization  error  dominates  NRBC  error  in  the  low 
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order  element  (order  1  and  2)  cases.  However,  by  increasing  the  polynomial  order  of  the 
elements,  this  spatial  discretization  error  decreases  and  allows  the  true  accuracy  of  the 
NRBC  to  be  realized. 


Figure  15:  Semi-infinite  channel  L'q  error  versus  NRBC  and  spectral  element  order.  Do¬ 
main  is  discretized  into  9,409  points  for  all  spectral  element  orders  with  zero  advection  at 

t  =  3. 


B.  SEMI-INFINITE  HORIZONTAL  CHANNEL 

For  most  problems,  where  initial  data  is  generated  from  physical  measurements, 
it  would  be  impossible  to  generate  an  analytic  solution  with  which  to  compare  the  NRBC 
solution.  In  this  experiment,  we  use  the  initial  data  as  described  in  (VI.  1)  and  (VI. 2)  to  gen¬ 
erate  the  waves  in  the  solution.  To  see  the  effect  of  the  boundary  condition,  we  compare 
our  solution  to  one  computed  on  a  larger  domain,  i.e.,  —  2  <  x  <  10  where  a  homogeneous 
Dirichlet  boundary  condition  /i(10,  y,  t)  —  0  is  prescribed  for  F E,  replacing  the  NRBC.  For 
this  experiment,  the  discretization  is  chosen  to  maintain  a  mesh  of  28,033  grid  points  for 
each  spectral  order.  Time  integration  is  performed  with  RK-8  to  ensure  the  time  discretiza¬ 
tion  is  not  a  limiting  factor  in  computing  the  reference  solution.  We  then  solve  the  extended 
domain  solution  for  t  =  3,  ensuring  that  the  disturbance  has  time  to  propagate  through  the 
artificial  boundary,  yet  has  not  had  time  to  reach  the  eastern  boundary. 
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1. 


Semi-Infinite  Channel  with  Zero  Advection  Results 


In  Figure  16(a),  we  plot  the  reference  solution  on  the  top  panel  and  the  solution  of 
the  truncated  domain  using  the  J  =  4  G-N  NRBC  on  the  bottom  panel.  The  solution  on  the 
truncated  domain  uses  4th  order  spectral  elements  on  a  24  x  24  element  grid,  discretizing 
the  domain  into  9,409  grid  points.  Qualitatively  speaking,  the  results  show  very  little  re¬ 
flection  using  the  combination  of  high-order  elements  and  NRBC.  Quantitative  results  can 


Expanded  Domain  Reference  Solution 


NRBC  Solution  (J=4) 


(a)  h(x,y,  3):  U  =  0,V  =  0 


(b)  £2  :  u  =  o,  V  =  0 


Figure  16:  Semi-infinite  channel,  4th  order  spectral  elements  (J  =  4)  using  cosine  pulse 
initial  condition  and  zero  advection.  Left  Plot:  Contour  plot  showing  h(x,  y,  3)  on  ex¬ 
tended  and  truncated  domains.  Right  Plot:  Corresponding  error  versus  NRBC  and 
spectral  element  order.  Domain  is  discretized  into  9,409  points  for  all  spectral  element 
orders. 


be  observed  in  Figure  16(b)  showing  the  error  on  f)  as  a  function  of  SE  and  NRBC  order 
(J  =  1, . . . ,  10, 15,  20).  The  number  of  elements  is  adjusted  for  each  polynomial  order  to 
maintain  an  equal  number  of  points  (9,409)  that  the  domain  is  discretized  into. 

It  is  clear  that  increasing  the  NRBC  order  yields  significant  gains  in  accuracy,  but 
by  J  =  5,  the  spatial  discretization  error  dominates  NRBC  error  in  the  low  order  element 
(order  1  and  2)  cases.  However,  by  increasing  the  spectral  order  of  the  elements,  this  spatial 
discretization  error  decreases  and  allows  the  true  accuracy  of  the  NRBC  to  be  realized. 
Similar  qualitative  and  quantitative  results  are  found  using  the  Gaussian  initial  data  and  are 
shown  in  Figure  17. 
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Expanded  Domain  Reference  Solution 


(a)  h(x,  y,  3):  U  =  0,  V  =  0  (b)  L2n:  U  =  0,  V  =  0 


Figure  17:  Semi-infinite  channel,  Ath  order  spectral  elements  ( J  =  4)  using  Gaussian  initial 
condition  with  zero  advection.  Left  Plot:  Contour  plot  showing  h(x,y,  3)  on  extended 
and  truncated  domains.  Right  Plot:  Corresponding  L‘f}  error  versus  NRBC  and  spectral 
element  order.  Domain  is  discretized  into  9,409  points  for  all  spectral  element  orders. 

2.  Semi-Infinite  Channel  with  Constant  Advection  Results 

In  this  section,  we  continue  the  analysis  on  the  semi-infinite  channel,  this  time 
adding  constant  advection  in  various  directions.  In  this  setting,  we  again  must  consider 
the  selection  of  the  Higdon  parameters  Cj.  For  the  semi-infinite  channel,  we  make  our 
educated  choice  described  in  Chapters  III.B  and  IV.D  of  augmenting  the  parameter  value 
with  the  added  advection.  In  the  case  of  the  eastern  boundary  NRBC,  this  yields  a  choice 
of  Cj  =  c0  +  U  which  has  the  simplifying  feature  of  making  the  boundary  condition  only 
first  order  in  time. 

In  Figure  18(a),  we  replicate  the  comparison  plot  between  the  reference  solution 
and  the  truncated  domain  solution.  Here,  we  use  J  =  4  G-N  NRBC  and  Ath  order  spectral 
elements  on  the  same  24  x  24  element  grid,  this  time  with  constant  advection  U  =  0.1 
and  V  =  0.  Qualitatively  speaking,  the  results  show  very  little  reflection  using  the  com¬ 
bination  of  high-order  elements  and  NRBC.  Further,  if  comparing  this  solution  with  the 
zero  advection  case,  one  notes  the  expansion  of  the  wave  in  the  direction  of  advection  and 
compression  where  the  wave  is  traveling  against  the  direction  of  advection.  These  results 
match  intuition  and  is  the  same  behavior  that  van  Joolen  noted  in  [16,  p.  108]. 
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Expanded  Domain  Reference  Solution 


(a)  h(x,y,  3):  U  =  0.1,  V  =  0 


(b)  L\{.  U  =  0.1,  V  =  0 


Figure  18:  Semi-infinite  channel,  4th  order  spectral  elements  (J  =  4)  using  cosine  pulse 
initial  condition  with  advection  velocities  U  =  0.1  and  V  =  0.  Left  Plot:  Contour  plot 
showing  h(x,  y,  3)  on  extended  and  truncated  domains.  Right  Plot:  Corresponding  L2Q 
error  versus  NRBC  and  spectral  element  order.  Domain  is  discretized  into  9,409  points  for 
all  spectral  element  orders. 


Quantitative  results  can  be  observed  in  Figure  18(b)  showing  the  error  on  0  as 
a  function  of  SE  and  NRBC  order  ( J  =  1, . . . ,  10, 15,  20).  The  number  of  elements  is 
again  adjusted  for  each  polynomial  order  to  maintain  an  equal  number  of  points  (9,409) 
that  the  domain  is  discretized  into.  It  is  again  clear  that  increasing  the  NRBC  order  yields 
significant  gains  in  accuracy,  but  the  spatial  discretization  error  dominates  NRBC  error  in 
the  low  order  element  (order  1  and  2)  cases. 

Error  is  monotonically  decreasing  for  each  SE  order,  however,  an  interesting  “dip” 
occurs  for  J  =  3.  Anomalies  such  as  this  can  occur  when  choosing  the  C/s  in  a  general 
manner  as  we  have  undertaken  in  this  example.  While  the  reflection  coefficient  guarantees 
that  the  reflection  caused  by  the  boundary  will  decrease  as  J  increases,  it  says  nothing  about 
how  much  it  will  decrease.  This  depends  heavily  on  the  choice  of  the  C/s;  in  this  example, 
it  is  believed  that  our  general  C3  choice  for  J  =  3  happened  to  annihilate  a  significant 
wave  mode  or  modes,  resulting  in  better  performance  of  the  boundary  condition.  Similar 
qualitative  and  quantitative  results  are  found  using  the  Gaussian  initial  data  and  are  shown 
in  Figure  19. 
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Expanded  Domain  Reference  Solution 


NRBC  Solution  (J=4) 


(a)  h(x,y,  3):  U  =  0.1,  V  =  0 


(b)  L2n:  U  =  0.1,  V  =  0 


Figure  19:  Semi-infinite  channel,  Ath  order  spectral  elements  ( J  =  4)  using  Gaussian  initial 
condition  with  advection  velocities  U  —  0.1  and  V  —  0.  Left  Plot:  Contour  plot  showing 
h(x,  y,  3)  on  extended  and  truncated  domains.  Right  Plot:  Corresponding  L'A  error  versus 
NRBC  and  spectral  element  order.  Domain  is  discretized  into  9,409  points  for  all  spectral 
element  orders. 

As  discussed  in  Chapter  III,  the  Higdon  NRBC  implicitly  assumes  that  any  wave 
impinging  on  the  NRBC  is  traveling  primarily  as  a  plane  wave  normal  to  the  boundary.  The 
previous  example,  where  the  advection  velocity  was  in  the  same  direction  as  the  channel 
does  not  significantly  challenge  this  assumption.  In  other  words,  to  really  test  the  value  of 
the  boundary  condition,  one  must  try  advection  velocities  with  some  tangential  component 
to  the  boundary.  For  these  experiments,  we  consider  only  the  Gaussian  initial  condition8. 

Figures  20(c)  and  20(d)  show  the  same  L‘iA  plots  for  advection  velocities  in  other 
directions,  one  being  the  contrived  case  where  the  advection  is  perfectly  tangential  to  the 
NRBC.  These  LA  plots  correspond  with  the  contour  plots  directly  above  that  show  the 
comparative  solutions.  Examining  these  results,  it  is  clear  that  the  boundary  condition  - 
even  when  put  to  a  test  with  a  wave  pulse  containing  a  significant  tangential  component 
to  the  boundary  -  performs  well.  It  is  noted  that  the  order  of  the  error  suffers  more  under 
diagonal  advection  when  compared  to  its  individual  axial  counterparts.  This  may  be  due  to 

8Since  the  cosine  pulse  initial  condition  was  constructed  to  ensure  that  the  boundary  condition  oil  I’  y  and 
r s  were  automatically  satisfied,  any  tangential  advection  would  cause  this  condition  to  break  down.  For  this 
reason,  we  choose  to  illustrate  general  results  using  only  the  Gaussian  initial  condition. 
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the  additional  terms  activated  in  the  interior  and  boundary  formulations  (IV.  9)  and  (IV.  10) 
when  U  and  V  are  simultaneously  nonzero. 

It  should  be  further  noted  that  as  advection  velocities  are  taken  larger,  the  associated 
errors  grow  as  well.  In  fact,  experimentation  shows  that  the  formulation  actually  becomes 
unstable  for  high-order  NRBCs  as  the  velocities  approach  the  reference  wave  speed  c0. 
This  is  the  same  behavior  noted  by  van  Joolen  in  [31]  when  examined  in  a  finite  difference 
setting.  This  behavior  does  not  seem  to  manifest  itself,  however,  unless  the  advection 
velocities  are  taken  much  larger  than  would  be  expected  in  a  geophysical  setting.  Further, 
this  behavior  seems  to  be  exacerbated  in  all  diagonal  advection  cases  when  using  high- 
order  (N  =  16)  basis  functions.  To  counter  this  high-order  instability,  we  implement  high- 
frequency  wave  filtering  as  described  by  Giraldo  et  al.  in  [47,  50]  when  using  N  =  16 
basis  functions.  As  in  [47],  we  apply  the  filter  every  10  time-steps  at  20%  strength. 

C.  INFINITE  HORIZONTAL  CHANNEL 


For  the  next  set  of  experiments,  the  domain  is  an  infinite  channel  with  the  NRBCs 
located  at  x  =  ±2.  The  set-up  is  similar  to  that  of  the  semi-infinite  channel  in  that  no¬ 
flux  boundary  conditions  are  specified  on  Tat  and  Y 5.  Further,  T ^  is  exactly  the  same 
NRBC  defined  previously.  This  time,  however,  we  prescribe  another  NRBC  B  on  Yw- 
This  specific  situation  is  shown  in  Figure  21. 

The  adjustments  required  for  the  western  NRBC  follow  from  the  original  derivation 
of  the  Higdon  boundary  condition  as  given  in  Chapter  III.A.  Now,  instead  of  considering 
the  right  moving  component  of  the  wave  approaching  Ye,  we  consider  the  left  moving 
component  of  the  wave  approaching  Yw.  To  perfectly  absorb  the  wave  impinging  on  the 
boundary,  we  insist  that  F  =  0  in  (III.  10)  such  that  the  boundary  satisfies  the  one-way 
advection  equation  h(x,t )  =  G^x  +  (c0  —  U)t^j.  The  corresponding  Higdon  boundary 
condition  is  then  given  by: 


H.j  : 


0  on  r  iy  . 


(VI. 6) 
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Expanded  Domain  Reference  Solution  Expanded  Domain  Reference  Solution 


(a)  h(x,  y,  3):  U  =  0.0,  V  =  0.1  (b)  h(x,  y,  3):  U  =  0.1,  V  =  0.1 


(c)  Lq:  U  =  0.0,  V  =  0.1 


(d)  L2n\  U  =  0.1,  V  =  0.1 


Figure  20:  Semi-infinite  channel,  4th  order  spectral  elements  ( J  =  4)  using  Gaussian  initial 
condition  with  advection  velocities  specified.  Top  Plots:  Contour  plots  showing  h(x,  y,  3) 
on  extended  and  truncated  domains.  Bottom  Plots:  Corresponding  Li,  error  versus  NRBC 
and  spectral  element  order.  Domain  is  discretized  into  9,409  points  for  all  spectral  element 
orders. 


In  the  experiments  that  follow,  we  use  the  initial  data  as  described  in  (VI.  1)  and 
(VI. 2)  to  generate  the  waves  in  the  solution.  To  see  the  effect  of  the  boundary  condi¬ 
tion,  we  compare  our  solution  to  one  computed  on  a  larger  domain,  i.e.,  —  6  <  x  <  6 
where  homogeneous  Dirichlet  boundary  conditions  h(— 6,  y,  t)  —  0  and  h( 6,  y,t)  =0  are 
prescribed  for  Tw  and  TE,  replacing  the  NRBCs.  The  discretization  is  again  chosen  to 
maintain  a  mesh  of  28,033  grid  points  for  each  SE  order.  Time  integration  is  performed 
with  RK-8  to  ensure  the  time  discretization  is  not  a  limiting  factor  in  computing  the  ref¬ 
erence  solution.  We  then  solve  the  extended  domain  solution  for  t  —  3,  ensuring  that  the 
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yN=2 


ys=-2 


Figure  21:  The  infinite  channel  domain  under  consideration.  Domain  is  truncated  by  artifi¬ 
cial  boundaries  at  xw  =  —2,  and  xe  =  2. 


disturbance  has  time  to  propagate  through  the  artificial  boundaries,  yet  has  not  had  time  to 
reach  the  Dirichlet  boundaries. 

1.  Weak  Form  Adjustments 

Using  a  similar  strategy  for  introducing  a  set  of  auxiliary  variables  for  the  western 
NRBC,  applying  them  to  the  KGE  equivalent,  then  converting  any  normal  derivatives  on 
the  boundary  to  time  and  tangential  boundaries,  results  in  another,  very  similar  formulation 
that  is  directly  incorporated  into  the  weak  form  (IV.3).  The  details  of  this  formulation  can 
be  viewed  in  Appendix  G.  The  selection  of  parameters  Cj  follows  the  “convenient”  choice 
as  developed  for  the  eastern  boundary,  namely,  to  remove  the  second  order  in  time  auxiliary 
variable  term  by  augmenting  the  reference  wave  speed  with  advection.  This  choice  is 
Cj  =  Co  —  U  for  the  western  boundary.  Similar  experiments  to  those  run  in  the  semi¬ 
infinite  channel  were  conducted  using  various  advective  velocities. 

2.  Infinite  Channel  with  Various  Advection  Velocities  Results 

Qualitative  and  quantitative  results  using  the  cosine  pulse  initial  condition  are  shown 
in  Figure  22.  In  Figure  22(a),  we  plot  the  reference  solution  on  the  top  panel  and  the  solu¬ 
tion  of  the  truncated  domain  using  the  J  =  4  G-N  NRBC  on  the  bottom  panel.  Quantitative 
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results  can  be  observed  in  Figure  22(b)  showing  the  error  on  as  a  function  of  SE  and 
NRBC  order.  In  Figure  22(c),  we  replicate  the  comparison  plot  between  the  reference  solu¬ 
tion  and  the  truncated  domain  solution,  this  time  with  left  to  right  advection.  Quantitative 
results  can  be  observed  in  Figure  18(b)  showing  the  error  on  Q  as  a  function  of  SE  and 
NRBC  order.  In  both  cases,  the  number  of  elements  is  again  adjusted  for  each  polynomial 
order  to  maintain  an  equal  number  of  points  (9,409)  that  the  domain  is  discretized  into. 


Expanded  Domain  Reference  Solution 


-3-2-101234 


NRBC  Solution  (J=4) 


(a)  h{x,y,  3):  U  =  0,V  =  0 


(b)  L2n:  U  =  0,  V  =  0 


Expanded  Domain  Reference  Solution 


|  • 


NRBC  Solution  (J=4) 


(c)  h(x,y,  3):  U  =  0.1,  V  =  0 


(d)  L2n:  U  =  0.1,  V  =  0 


Figure  22:  Infinite  channel,  4th  order  spectral  elements  (J  =  4)  using  cosine  pulse  initial 
condition  with  advection  velocities  specified.  Left  Plots:  Contour  plots  showing  h(x,  y,  3) 
on  extended  and  truncated  domains.  Right  Plots:  Corresponding  L'i}  error  versus  NRBC 
and  spectral  element  order.  Domain  is  discretized  into  9,409  points  for  all  spectral  element 
orders. 
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Similar  qualitative  and  quantitative  results  are  found  using  the  Gaussian  initial  data 
and  are  shown  in  Figure  23.  It  is  noted  that  the  order  of  the  error  again  suffers  more  under 
diagonal  advection  when  compared  to  its  individual  axial  counterparts. 

D.  OPEN  DOMAIN  CONSIDERATIONS 

In  the  construction  of  the  G-N  NRBCs  presented  thus  far,  we  have  assumed  that  all 
boundaries  are  aligned  with  the  axial  coordinates.  Further,  no  two  NRBCs  have  ever  been 
placed  adjacent  to  each  other.  In  this  section,  we  examine  the  consequences  of  placing 
NRBCs  adjacent  to  each  other  in  two  configurations,  namely  the  quarter  plane  and  the 
open  plane.  This  set-up  is  similar  to  the  channel  configurations  described  thus  far  in  that 
the  infinite  domain  is  truncated  via  artificial  boundaries  B  thus  dividing  the  domain  into  a 
finite  computational  domain  f)  and  a  residual  domain  D.  The  only  thing  that  changes  is  the 
configuration  of  the  artificial  boundaries. 

Specifically,  the  quarter  plane  is  described  as  a  domain  that  is  bounded  by  physical 
boundaries  Tw  and  Ts.  NRBCs  are  introduced  at  x  =  xE  and  y  =  yN.  The  physical 
boundaries  are  homogeneous  Dirichlet  conditions  h  =  0  on  and  Ts.  This  set-up  is 
illustrated  in  Figure  24(a).  The  open  plane  is  described  as  a  domain  that  is  unbounded  on 
all  sides.  To  compute  a  solution  on  such  a  domain,  NRBCs  are  introduced  at  x  =  %,  xE 
and  y  =  ys,  yN-  This  setup  is  illustrated  in  Figure  24(b).  Artificial  boundaries  for  F^  and 
T N  are  developed  as  outlined  in  Appendix  G. 

1.  Corner  Compatibility  Concerns 

To  begin  this  discussion,  consider  the  quarter  plane.  A  source  of  concern  is  the 
method  of  handling  the  intersection  point  of  VE  and  T  After  all,  the  auxiliary  variable 
form  described  in  (III. 20)  and  (G.9)  are  PDEs  themselves  and  therefore  require  appropriate 
boundary  conditions  to  be  well-posed.  In  the  channel,  the  no-flux  conditions  specified  by 
the  problem  statement  were  applied  to  the  auxiliary  variables  to  make  the  problem  well- 
posed.  In  the  quarter  plane,  we  have  the  homogeneous  Dirichlet  conditions  for  the  western 
boundary  of  Tat  and  the  south  boundary  of  T E,  but  there  are  no  such  boundary  conditions 
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Expanded  Domain  Reference  Solution 


K 


□ 


NRBC  Solution  (J=4) 


(a)  h(x,y,  3):  U  =  0,V  =  0 


Expanded  Domain  Reference  Solution 


■ 


NRBC  Solution  (1=4) 


(c)  h(x,y,  3):  U  =  0.1,  V  =  0 


Expanded  Domain  Reference  Solution 


tc — =■>! 


NRBC  Solution  (J=4) 


(e)  h(x,y,  3):  U  =  0,V  =  0.1 


Expanded  Domain  Reference  Solution 


(g)  h(x,y,3):  U  =  0.1,  V  =  0.1 


(b)  Lfl:U  =  0,  V  =  0 


(d)  L2n:  U  =  0.1,  V  =  0 


(f)  Lbu  =  0,V  =  0.1 


(h)  L2n:  U  =  0.1,  V  =  0.1 


Figure  23:  Infinite  channel,  4//l  order  spectral  elements  (J  =  4)  using  Gaussian  initial 
condition  with  advection  velocities  specified.  Left  Plots:  Contour  plots  showing  h(x,  y,  3) 
on  extended  and  truncated  domains.  Right  Plots:  Corresponding  L q  error  versus  NRBC 
and  spectral  element  order.  Domain  is  discretized  into  9,409  points  for  all  spectral  element 
orders.  „„ 


(a)  Quarter  Plane  Domain 


(b)  Open  Plane  Domain 


Figure  24:  Left  Plot:  A  semi-infinite  domain  Q  truncated  by  artificial  boundaries  T E  and 
T at.  Right  Plot:  An  infinite  domain  Q  truncated  by  artificial  boundaries  Us,  Fr-  T  v  and 

r  w- 

that  can  be  applied  to  the  east  of  Tn  and  north  of  Y E.  The  question  therefore  arises,  what 
are  the  appropriate  boundary  conditions  for  the  auxiliary  variables  at  these  points?  This 
problem  is  compounded  in  the  open  domain  as  there  are  no  explicitly  defined  boundary 
conditions  for  any  of  the  boundaries  of  the  auxiliary  variables. 

2.  Use  of  Sommerfeld  Radiation  Boundary  Conditions  for  Auxiliary  Vari¬ 
able  Boundary  Conditions 

For  this  analysis,  we  suggest  that  the  desired  behavior  of  the  boundary  data  on 
the  auxiliary  variables  at  these  comers  should  minimize  auxiliary  variable  reflection  back 
into  the  computational  (boundary)  domain.  In  other  words,  the  auxiliary  variables  should 
themselves  be  non-reflecting.  Ultimately,  we  would  like  these  boundary  conditions  to  be 
easily  implementable  while  still  maintaining  the  true  essence  of  the  auxiliary  variables. 
To  implement  this  behavior ,  we  consider  a  simple  order  J  =  1  G-N  NRBC  (Sommerfeld 
condition)  for  the  intersection  points  of  two  NRBCs,  i.e., 

4>'AxE,yN)  =  — —  <j>j(xE,yN)  forrB  (VI.7) 

c0,y 

1  . 

<j>j{xE,yN)  = - 4>j{xE,yN)  for  Tat  (VI.8) 

Q),# 
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Here,  c0)J/  and  c0tX  are  the  y  and  x  components  of  the  radial  wave  velocity  and  prime  indi¬ 
cates  the  tangential  derivative  along  the  particular  boundary,  i.e.,  0'-  =  0^-  along  T E  and 
0'-  =  along  T at.  We  now  recall  the  auxiliary  form  of  the  boundary  condition  for  the  0.; 

terms  onfE 


aj  /  Ci0,'-1  dTE+Kj  /  Ciftj-l  dF E  —  Xy  /  Ciftj-l  dF E  +  Pj  /  C i4>jdTE 

Jye  Jte  drE  drE 

-7  [  Wj  drE  -f  f  Q&j-i  dr E  =  \x  f  Ci^j+i  drE- 

JvE  JrE  drE 

To  ensure  that  the  auxiliary  form  lies  in  H 1  (T E),  we  integrate  the  second  order  in  space 
term  by  parts  to  yield 


dr E+Hj  drE  +  \  /  Cti-i  drE  -  XyQti 


Vn 


+  Pj  /  00,  dr 7 


2/S 


-7  /  00',- drE-/2  /  Q0j-i  drB  =  Ax  /  c*0j+i  rfTi 


We  now  see  that  the  boundary  term  contains  0'  evaluated  at  the  northern  boundary.  To 
implement  the  non-reflecting  behavior,  we  make  the  substitution  (VI.7)  into  the  boundary 
term.  The  complete  weak  boundary  form  (for  T E)  then  takes  the  form 


ay  /  00,-1  dr  e+Kj  [  drE  +  Xy  [  00 7-1  drE  +  Xy—Q^j-i  +  Pj  [  Ci4>j  drE 

JrE  J  rE  JrE  co,y  JrE 

— 7  [  0 4>'j  dr E  —  f2  f  00, -i  dr E  =  xx  f  00,+ 1  dr E 

JrE  JrE  drE 

for  j  =  1, . . . ,  J  —  1;  V  0,  0,'  €  VrB  and  <pj  =  0  at  ?/s.  A  similar  construct  is  readily 
computed  for  T at. 

The  only  thing  left  to  do  is  consider  the  problem  of  “double  counting”  the  contri¬ 
bution  at  the  comer.  In  other  words,  there  are  two  values  for  the  auxiliary  variables  at  the 
comer;  the  one  resulting  from  the  evaluation  of  T ^  and  the  one  from  T N-  Which  aux¬ 
iliary  variable  contribution  at  the  corner  should  be  used,  that  of  r E  or  that  of  T  Y?  For 
this  analysis,  we  adopt  the  “node-splitting”  approach  described  by  Pozrikidis  [44,  p.  215], 
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which  amounts  to  averaging  the  corner  auxiliary  variable  values.  Of  course,  this  formula¬ 
tion  is  inherently  low  order  in  the  boundary  treatment  of  the  auxiliary  variables.  As  such, 
we  would  not  expect  to  see  spectral  convergence  as  shown  in  previous  examples,  however, 
there  should  be  improvement  over  the  J  =  1  Sommerfeld  condition. 

3.  Corner  and  Open  Domain  with  Zero  Advection  Results 

Figure  25  shows  a  series  of  contour  plots  showing  how  the  initial  disturbance  prop¬ 
agates  through  the  domain  for  0  <  t  <  4.5  with  zero  advection.  In  this  example,  we  run  the 
simulation  using  4th  order  basis  functions  on  a  24  x  24  -  element  grid  (9,409  global  points), 
using  NRBC  order  J  =  4.  The  simulation  is  run  just  long  enough  for  the  primary  wave 
to  exit  the  computational  domain.  Qualitatively  speaking,  the  results  appear  to  behave  as 
desired  -  allowing  the  wave  to  propagate  through  the  NRBC  unimpeded. 

Time  Evolution  of  the  Quarter  Plane  Gaussian 


t=0.25  t=0.5  t=0.75  t=l  t=l  .25  t=1.5 


Figure  25:  Time  Evolution  of  quarter  plane  Gaussian  (NRBC  on  and  T E)  using  4th 

order  spectral  elements  ( J  =  4)  with  zero  advection. 
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Quantitatively,  the  results  confirm  that  errors  are  not  exponentially  decaying  as  a 
function  of  NRBC  order.  As  the  order  of  the  NRBC  is  increased,  the  crude  approximation 
of  the  corner  boundary  condition  on  the  auxiliary  variables  shows  its  weakness.  In  fact, 
experimentation  shows  that  the  boundary  condition  error  quickly  overtakes  spatial  and  time 
discretization  errors.  Taking  the  Gaussian  initial  condition  with  tf  —  3.0  for  various  basis 
function  orders  and  boundary  condition  orders  yields  the  L\  errors  (using  an  extended 
domain  solution  as  the  reference)  as  shown  in  Table  2. 

Table  2:  Lq  Error  as  a  function  of  NRBC  Order  for  quarter  plane  using  various  spectral 
element  orders.  Gaussian  initial  condition  and  zero  advection  is  used. 


NRBC  Order 

Error 

Finear  Elements 

L  f)  Error 
Order  2  Elements 

L'o  Error 
Order  4  Elements 

J  =  1 

0.11002 

0.10519 

0.10470 

J  =  2 

0.07204 

0.04173 

0.05126 

J  =  3 

0.04067 

0.03377 

0.03277 

J  =  4 

0.03948 

0.02919 

0.03274 

J  =  5 

0.03647 

0.02579 

0.02941 

J  =  10 

0.03446 

0.02517 

0.02539 

J  =  20 

0.03446 

0.02513 

0.02526 

Similar  experiments  were  performed  for  the  open  plane.  G-N  boundary  conditions 
are  implemented  along  all  four  boundaries  and  the  order  1  Sommerfeld  boundary  condition 
is  applied  to  each  auxiliary  variable  boundary  as  described  in  the  previous  section.  Figure 
26  shows  a  series  of  contour  plots  showing  how  the  initial  disturbance  propagates  through 
the  domain  for  0  <  t  <  4.5  with  zero  advection.  In  this  example,  we  run  the  simulation 
using  4th  order  basis  functions  on  a  24  x  24  -  element  grid  (9,409  Global  Points),  using 
NRBC  order  J  =  4.  Again,  qualitatively  speaking,  the  results  appear  to  behave  as  desired 
-  allowing  the  wave  to  propagate  through  the  NRBC  unimpeded. 

Again,  taking  the  Gaussian  initial  condition  with  tf  —  3.0  for  various  basis  function 
and  boundary  condition  orders  yields  the  L\  errors  (using  an  extended  domain  solution  as 
the  reference)  as  shown  in  Table  3.  Two  main  observations  can  be  drawn  from  these  results. 
First,  the  major  source  of  error  appears  to  be  with  the  boundary  treatment.  While  there  are 
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Time  Evolution  of  the  Open  Plane  Gaussian 


t=0.25  t=0.5  t=0.75  t=l  t=  1 .25  t=1.5 


Figure  26:  Time  Evolution  of  open  plane  Gaussian  (NRBC  on  all  boundaries)  using  4th 
order  spectral  elements  ( J  =  4)  with  zero  advection. 

modest  gains  made  by  increasing  the  spectral  element  order,  once  J  =  3,  the  errors  for 
all  spectral  element  orders  are  nearly  the  same.  Second,  even  though  the  Liy  results  are 
less  impressive  than  the  channel  experiments  for  high-order  NRBCs,  there  is  significant 
improvement  from  J  =  1  (Sommerfeld  condition)  to  higher  J  -  even  if  the  improvement 
is  far  from  exponential. 

With  this  said,  one  may  be  concerned  with  the  stability  and  long  term  behavior 
of  this  NRBC  scheme  employed  in  the  quarter  and  open  domain  planes.  Of  course,  as 
t  — >■  oo,  one  would  expect  h  0.  To  gain  a  quantitative  handle  on  this,  consider  the 
oo-norm  defined  as  follows: 

IN  loo  =  max  \hi\ 
l<i<  Np 

where  Np  is  the  number  of  points  in  Q.  We  choose  this  norm  simply  to  get  an  estimate  of 
how  much  of  the  initial  disturbance  is  left  in  the  computational  domain  after  a  substantial 
amount  of  time  has  passed.  Using  our  now  standard  test  case  of  4th  order  spectral  elements 
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Table  3:  L ^  Error  as  a  function  of  NRBC  Order  for  open  plane  using  various  spectral 
element  orders.  Gaussian  initial  condition  and  zero  advection  is  used. 


NRBC  Order 

L'ij  Error 
Linear  Elements 

L\  Error 
Order  2  Elements 

L’o  Error 
Order  4  Elements 

J  =  1 

0.23695 

0.21538 

0.21102 

J  =  2 

0.09145 

0.08376 

0.04345 

J  =  3 

0.06494 

0.04056 

0.04144 

J  =  4 

0.06466 

0.03880 

0.03837 

J  =  5 

0.06454 

0.03873 

0.03776 

J  =  10 

0.06441 

0.03794 

0.03767 

J  =  20 

0.06441 

0.03794 

0.03767 

on  a  24  x  24  element  grid  with  NRBC  order  J  —  4,  when  computed  for  t  =  1000  with 
J  =  10,  it  was  found  that  |  |/i|  |oo  =  1-03  x  10~17forthe  quarter  plane  and  |  |/i|  |oo  =  5.30  x 
10~19  for  the  open  domain;  in  both  cases,  essentially  zero  throughout  the  computational 
domain.  While  this  is  not  a  rigorous  stability  analysis,  it  does  experimentally  suggest  a 
stable  formulation. 

4.  Corner  and  Open  Plane  Domain  with  Constant  Advection  Results 

It  can  be  shown  (see  Appendix  H  for  details)  that  when  working  in  the  open  plane,  to 
examine  the  behavior  of  the  KGE  under  constant  advective  velocities  in  various  directions, 
a  simple  change  of  coordinate  system  can  recast  the  problem  into  a  much  simpler  problem. 
The  simplified  problem  is  one  where  advection  is  in  only  one  direction  aligned  with  the 
new  coordinate  system.  This  implies  that  when  examining  the  open  plane  under  advection, 
it  is  sufficient  to  test  only  cases  where  advection  is  in  the  x  or  y  direction.  The  benefits  of 
this  change  of  coordinate  system  include  reducing  the  computational  overhead,  as  well  as 
minimizing  various  errors  due  to  the  more  complex  formulation  if  viewed  in  the  original 
coordinate  system.  It  should  be  noted,  however,  that  the  formulation  discussed  in  VI.D.3 
still  results  in  a  stable  formulation  when  diagonal  advection  is  applied  to  the  solution.  To 
see  this,  examine  Figure  27  where  north-east  advection  is  applied  and  the  co,y  and  c0,;,:  terms 
have  been  adjusted  by  advection. 
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Time  Evolution  of  the  Quarter  Plane  Gaussian 


t=0.25  t=0.5  t=0.75  t=l  t=l  .25  t=1.5 


Figure  27:  Time  Evolution  of  open  plane  Gaussian  (NRBC  on  all  boundaries)  using  4th 
order  spectral  elements  (J  =  4)  with  advection  velocities  U  —  0.1  and  V  =  0.1. 


The  corresponding  L ^  errors  are  presented  for  various  SE  and  NRBC  orders  in 
Table  4. 

E.  EFFECTS  OF  HIGH-ORDER  TIME  INTEGRATION 

At  the  outset  of  this  work,  it  was  believed  that  at  some  point  the  improvements  real¬ 
ized  by  increasing  the  spatial  discretization  and  the  order  of  the  NRBC  would  eventually  be 
limited  by  the  time  integration  scheme  [52].  To  this  end,  the  order  of  the  time  integration 
scheme  was  varied  to  examine  the  effects  of  time  integration  on  accuracy  of  the  solution. 
As  has  already  been  presented,  gains  made  by  increasing  the  order  of  the  NRBC  halt  for 
lower  order  spectral  elements  after  J  =  5.  Early  experiments  showed  that  even  for  high 
order  (order  8  and  16)  spectral  elements,  the  gains  made  by  increasing  the  order  of  the 
NRBC  are  limited  at  some  point  using  classical  RK-4  time  integration. 
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Table  4:  Error  as  a  function  of  NRBC  Order  for  open  plane  using  various  spectral 

element  orders  .  Gaussian  initial  condition  with  advection  velocities  U  —  0.1  and  V  =  0.1 
used. 


NRBC  Order 

L  j,  Error 
Linear  Elements 

L\  Error 
Order  2  Elements 

L’o  Error 
Order  4  Elements 

J  =  1 

0.26914 

0.26529 

0.16998 

J  =  2 

0.09571 

0.08642 

0.03810 

J  =  3 

0.04391 

0.02711 

0.02631 

J  =  4 

0.04061 

0.02514 

0.02467 

J  =  5 

0.04028 

0.02413 

0.02460 

J  =  10 

0.04021 

0.02410 

0.02460 

J  =  20 

0.04018 

0.02409 

0.02460 

For  this  experiment,  we  consider  the  KGE  on  a  semi-infinite  channel  with  hFw  =  0. 
To  ensure  that  any  boundary  or  time  effects  are  not  masked  by  the  interior  discretization, 
24th  order  spectral  elements  are  used  on  a  fine  mesh  consisting  of  4,753  global  points.  The 
Gaussian  initial  condition  is  used  and  is  evaluated  until  t  —  4.  The  reference  solution  in 
this  case  was  computed  as  described  previously,  except  this  time  using  24th  order  spectral 
elements  on  a  fine  mesh  consisting  of  9,457  global  points.  Time  integration  was  performed 
with  a  10/,!  order  Runge-Kutta  scheme  using  a  time-step  chosen  to  ensure  a  Courant  number 
of  0.1. 

As  can  be  observed  in  Figure  28,  gains  made  by  improving  the  time  integration 
matter  only  if  combined  with  high-order  treatment  of  the  boundary.  Conversely  -  gains 
using  high-order  treatment  of  the  boundary  can  only  be  realized  if  there  is  a  high-order 
treatment  of  the  time  integration.  It  should  be  noted  that  these  results  (error  on  the  order  of 
10~10)  cannot  be  observed  unless  high-order  treatment  of  the  interior  also  accompanies  the 
high-order  treatment  of  the  boundary  and  time.  Several  experiments  were  conducted  that 
varied  the  order  of  the  interior,  boundary  and  time  integration  [41,  53].  The  clear  result 
was  that  without  high-order  treatment  of  all  components  in  concert,  convergence  to  the 
reference  solution  is  stalled. 
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Figure  28:  Error  in  the  SE  Solution  of  KGE  using  NRBCs  of  various  order  as  a  function  of 
time  integration  order. 

We  believe  as  a  practical  matter  that  the  order  of  all  components  of  the  numerical 
solution  (spatial,  boundary  and  time)  should  be  chosen  to  ensure  that  the  numerical  method 
is  the  strongest  (in  accuracy)  component  of  the  model.  If  high-order  treatment  of  any 
of  the  three  components  is  missing,  the  high-order  treatment  of  the  other  components  is 
essentially  wasted.  For  models  where  parameters  and  data  have  associated  measurement 
and  parameter  errors,  the  numerical  method  should  be  chosen  to  maintain  at  least  the  same 
accuracy. 
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VII.  TOWARDS  ARBITRARY  DOMAINS 


In  the  previous  chapters,  we  examined  the  G-N  boundary  formulation  using  spectral 
elements  on  unstructured  quadrilaterals.  The  physical  boundaries,  however  were  perfectly 
aligned  with  the  coordinate  system  axes.  This  was  convenient  since  it  allowed  the  use  of 
the  G-N  auxiliary  formulation.  The  power  of  spectral  elements  lies  not  only  in  its  ability  to 
compute  high-order  accurate  solutions,  but  also  in  its  ability  to  handle  complex  geometries. 
While  we  demonstrated  exponential  error  convergence  in  channel  experiments  when  high- 
order  treatment  was  applied  to  spatial,  boundary  and  time  components  of  the  problem,  this 
exponential  convergence  broke  down  when  applied  to  boundary  configurations  where  two 
NRBCs  were  adjacent  to  each  other.  In  short,  since  there  is  a  discontinuity  in  the  normal  at 
the  intersection  of  adjacent  NRBCs,  the  corner  was  the  problem. 

This  chapter  considers  what  happens  when  we  completely  remove  any  corners.  We 
first  examine  an  arbitrary  domain  where  the  boundary  can  be  of  any  shape.  It  will  be 
shown  that  there  are  insurmountable  complications  that  arise  when  using  the  G-N  auxiliary 
formulation  in  this  context.  In  this  case,  results  are  presented  for  various  domains  using  a 
first  order  non-reflecting  boundary  condition  and  high-order  G-N  when  certain  simplifying 
assumptions  are  made.  The  chapter  concludes  by  revisiting  a  boundary  condition  originally 
devised  in  1998  for  the  wave  equation  by  Hagstrom  and  Hariharan  and  modified  in  2003 
by  vanJoolen  et  al.  to  include  dispersion. 

A.  ARBITRARILY  SHAPED  BOUNDARIES 

Ideally,  we  would  like  to  directly  extend  the  work  presented  thus  far  to  remove 
the  problematic  “corner”  configuration  and  replace  it  with  a  continuous,  smooth,  closed 
boundary.  If  this  were  possible,  then  a  single  formulation  (instead  of  four  formulations 
combined  in  the  open  domain)  for  the  boundary  would  result.  The  benefits  to  this  type 
of  boundary  would  be  that  the  domain  could  be  “fit”  to  the  area  of  interest,  reducing  the 
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total  number  of  grid  points  required.  A  general  configuration  that  demonstrates  this  idea  is 
shown  in  Figure  29. 


Figure  29:  A  general  domain  0  truncated  by  artificial  boundary  F 


To  begin,  we  again  examine  the  Higdon  boundary  condition  of  order  J  as  first  pre¬ 
sented  in  Chapter  III  and  the  KGE  as  presented  in  Chapter  II  simplified  by  the  assumption 
of  zero  advection. 


Hj: 


h 


h  -  CqV2/?.  +  f2h 


0 

0 


on  r 


We  note  that  the  boundary  condition  and  the  KGE  are  described  in  two  different  coordinate 
systems  -  namely  (n,  r)  and  (x,  y )  respectively  where  n  and  r  are  the  normal  and  tangential 
directions  on  the  boundary.  If  we  consider  an  arbitrary  part  of  the  boundary  (T)  as  shown 
in  Figure  30,  we  can  find  a  way  to  express  the  standard  Cartesian  derivatives  in  terms  of 
normal  and  tangential  derivatives. 

Of  course,  in  the  most  general  case,  the  normal  and  tangential  vectors  are  depen¬ 
dent  on  the  position  on  the  boundary,  i.e.,  n  =  n{x,y)  and  t  =  ?{x,y).  These  normal 
components  can  be  computed  (see  Appendix  I)  for  a  particular  domain  by  considering  a 
change  of  coordinates  from  ( x ,  y)  to  (n,  r)  as  defined  by  the  linear  transformation  and  its 


82 


n 


Figure  30:  Components  of  normal  and  tangential  derivatives 
associated  differentiation  operator: 
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Since  the  transformation  is  necessarily  non-singular,  this  can  also  be  written  as: 
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(VII.  1) 


(VII.2) 


Now,  if  we  expand  the  Higdon  boundary  condition,  for  J  =  1,  we  get 


dnh  +  —h  =  0  =>  n  ■  V/i  =  —  h. 

Li  Li 


This  is  convenient  since  this  boundary  condition  can  be  directly  applied  to  the  KGE  zero 
advection  weak  integral  form 


/  Vih  dQ-c20  /  %n  ■  Vh  dT  +  cl  /  W;  ■  Vh  dtt  +  f  /  %h  dQ  =  0 

In  Jr  Jn  Jn 
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to  yield  the  first  order  Sommerfeld  formulation 


%h  d<l  - 

Li 


^ihdr  +  cl  /  Vtf*  ■  V/i  dfi  + /2  /  %hdn  =  0. 


(VII.  3) 


/n 


1.  Second  Order  (and  Higher)  Higdon  Boundary  Condition 

Recall  that  the  Higdon  boundary  condition  is  very  general.  It  can  be  applied  to  a 
variety  of  wave-type  problems  and  reflection  is  guaranteed  to  decrease  by  simply  increasing 
the  order  J.  They  suffer,  however,  from  an  implementation  point  of  view  since  there  are 
increasingly  high-order  spatial  and  temporal  derivatives.  Consider  a  second  order  Higdon 
boundary  condition 

#2  :  (dn  +  (dnh  +  -^h^j  =0  on  T. 

When  expanded,  this  boundary  condition  takes  the  form 

/  1  1  \  1  .. 
dri  n  J1  +  (  J  ^ n  ^  C  C  1  =  ^ 

Continuing  with  the  expansion  to  express  the  boundary  condition  in  terms  of  the  physical 
coordinate  system,  we  find  that  dnnh  is 

r\  r\ 

dnnh=—(dnh)  =  —(n-Vh) 
d  ^ 

~  {jixdxh  H-  rtydyh^j  n  •  V  (nxOxh  -j-  Tiydyhdj .  (VII. 4) 

The  key  point  to  take  away  from  (VII.4)  is  that  the  components  of  the  normal  vectors  are 
themselves  functions  of  x  and  y.  The  product  rule  dictates  that  we  must  then  compute 
the  x  and  y  derivatives  of  the  normals  in  order  to  yield  an  “exact”  representation  of  the 
higher  order  Higdon  boundary  condition.  This  expansion  has  two  direct  consequences  that 
undermine  efficient  implementation  of  the  formulation. 


84 


The  first  is  something  that  has  already  been  addressed,  namely  that  the  ability  to 
implement  the  model  is  challenged  -  especially  when  dealing  with  large  J .  For  each  order 
that  the  Higdon  formulation  increases,  high-order  derivatives  appear  for  h  as  well  as  the 
components  of  the  normal  vector.  The  second  consequence  is  even  more  problematic  in 
that  the  means  to  relate  the  boundary  formulation  back  into  the  interior  formulation  has 
been  lost.  In  the  case  of  the  Sommerfeld  condition  presented,  the  boundary  condition 
and  the  boundary  integral  term  (following  integration  by  parts  of  the  Laplacian  operator) 
matched  perfectly  -  thus  allowing  direct  substitution  of  the  boundary  condition  into  the 
interior  formulation. 

2.  G-N  on  the  Unstructured  Boundary 

The  G-N  formulation  was  designed  to  remedy  this  problem  of  increasingly  high- 
order  derivatives  by  recasting  them  into  a  system  of  low  order  derivatives.  If  we  try  this 
with  the  unstructured  boundary  formulation,  a  similar  auxiliary  form  to  that  presented  in 
Chapter  III.C)  is  obtained 

(dn  +  =0,-  j  —  (VII. 5) 

where  0O  =  h  and  0j  =  0.  Again,  the  function  h  satisfies  the  KGE  and  0i  is  obtained  by 
applying  the  linear  (constant  coefficient)  operator  (dn  +  to  h. 

Knowing  that  in  the  end,  we  would  like  to  have  an  auxiliary  variable  formulation 
that  contains  only  tangential  and  time  derivatives  (so  that  the  boundary  formulation  can  be 
discretized  only  on  the  boundary),  we  must  consider  the  equation  that  0,  satisfies.  It  can  be 
shown  that  when  the  KGE  on  the  boundary  is  recast  in  terms  of  the  normal  and  tangential 
coordinate  system  that  it  becomes  a  variable  coefficient  differential  equation  due  to  the 
presence  of  the  normal  components  nx(x,y),ny(x,y).  The  result  of  this  is  that  there  is 
no  general  KGE- like  equation  that  all  0,  ’s  will  satisfy.  In  fact,  every  time  we  increase  the 
order  of  the  boundary  condition,  additional  terms  such  as  those  encountered  in  (VII.4)  are 
accumulated.  In  the  end,  a  separate  formulation  that  contains  high-order  derivatives  will 
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have  to  be  devised  for  each  auxiliary  variable  introduced.  In  short,  this  implies  that  the 
G-N  auxiliary  variable  formulation  in  its  current  form  is  incompatible  with  an  unstructured 
boundary. 

To  handle  this  incompatibility,  we  consider  a  simplifying  assumption  that  the  curva¬ 
ture  on  the  boundary  is  small.  This  replicates  the  case  where  in  a  local  region,  the  boundary 
“looks”  like  a  straight  line  to  the  numerical  solution.  This  assumption  allows  for  all  deriva¬ 
tives  of  normal  components  to  be  neglected,  i.e., 


d_ 

dx 


ij^x) 


d  .  d  . 

^(»,)  =  ^K)  =  o. 


The  details  of  this  boundary  formulation  and  how  it  is  integrated  into  the  interior  scheme 
are  discussed  in  Appendix  I. 


B.  RESULTS  FOR  ADJUSTED  G-N  NRBCS  ON  ARBITRARY  DOMAINS 


Some  experiments  were  performed  using  the  (what  turned  out  to  be)  convergence 
limiting  assumption  of  small  curvature.  While  experiments  showed  stable  behavior  for  the 
zero  advection  case,  even  over  long  term  time  integrations,  the  convergence  was  again,  far 
from  exponential  in  nature.  For  the  first  set  of  experiments,  we  consider  how  the  formula¬ 
tion  outlined  in  (VII.3)  performs  for  various  boundary  shapes.  Admittedly,  this  formulation 
is  only  a  first  order  boundary  condition,  however  as  has  already  been  shown,  a  first  order 
condition  is  very  easy  to  implement  and  has  very  modest  computational  overhead.  The 
next  set  of  experiments  considers  the  adjusted  G-N  formulation.  For  these  experiments, 
we  consider  rectangular,  circular,  and  rounded  rectangular  boundaries  where  the  Gaussian 
initial  condition  is  used  to  generate  the  propagating  waves. 

Figure  3 1  shows  a  series  of  contour  plots  for  the  zero  advection  case  using  4lh  order 
spectral  elements  with  J  =  1  and  J  =  4.  Model  parameters  are  set  to  the  standard  test 
case  where  —  1,  f2  —  0.5  and  initial  data  as  described  in  (VI.2)  is  used  to  generate  the 
waves  in  the  solution.  To  see  the  effect  of  the  boundary  condition,  we  compare  our  solution 
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to  one  computed  on  a  larger  domain,  i.e.,  x,y  e  [—4,4]  where  a  homogeneous  Dirichlet 
boundary  condition  h  =  0  is  prescribed  for  T,  replacing  the  NRBC.  For  this  experiment, 
the  discretization  is  chosen  to  maintain  a  mesh  of  28,033  grid  points  for  each  polynomial 
order.  Time  integration  is  performed  with  RK-8  to  ensure  the  time  discretization  is  not  a 
limiting  factor  in  computing  the  reference  solution.  We  then  solve  the  extended  domain 
solution  for  t  —  3,  ensuring  that  the  disturbance  has  time  to  propagate  through  the  artificial 
boundary,  yet  has  not  had  time  to  reach  the  boundary.  The  number  of  elements  in  each  of 
the  NRBC  solutions  is  adjusted  to  ensure  approximately  3,  000  global  points  were  used. 

What  is  clear  from  these  plots  is  that  there  are  trade-offs  between  accurately  rep¬ 
resenting  the  G-N  NRBC  and  removing  the  problematic  corners.  Specifically,  the  square 
domain  does  not  have  to  make  an  approximation  for  small  curvature.  In  fact,  with  the 
exception  of  only  4  points  (comers)  in  the  global  domain,  the  G-N  NRBC  is  perfectly  rep¬ 
resented  by  the  arbitrary  domain  formulation.  This  is  why  there  is  significant  improvement 
between  the  J  —  1  and  J  =  4  cases.  The  rounded  square  and  circular  domains  see  less 
dramatic  improvement  as  J  is  increased.  While  the  problematic  comers  are  removed  in 
these  cases,  the  small  curvature  assumption  induces  error  in  the  G-N  NRBC,  thus  imposing 
a  convergence  bound  that  cannot  be  overcome  by  simply  increasing  the  order  of  the  NRBC. 

A  series  of  experiments  was  conducted  to  examine  the  normalized  L ^  error  as  a 
function  of  spectral  element  and  NRBC  order  for  each  of  the  NRBC  boundary  configu¬ 
rations.  What  is  clear  from  the  results  shown  in  Tables  5-7  is  that  the  errors  are  more  a 
function  of  the  NRBC  than  of  the  spectral  element  order.  In  short,  there  is  almost  no  gain 
observed  by  using  high-order  spectral  elements  since  much  of  the  error  is  generated  by  the 
NRBC.  When  compared  with  the  results  shown  in  Chapter  VI.D.3,  which  used  the  Som- 
merfeld  approximation  for  the  boundary  condition  of  the  auxiliary  variables,  the  results  are 
unimproved. 
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(a)  Ref.  Solution  (b)  Truncated  Ref.  Solution 


(c)  NRBC  J  =  1 


(d)  NRBC  J  =  1 


(e)  NRBC  J  =  1 


(f)  NRBC  J  =  4 


(g)  NRBC  J  =  4  (h)  NRBC  J  =  4 


Figure  31:  Open  Domain,  Ath  order  spectral  elements  using  Gaussian  initial  condition  with 
zero  advection.  Top  Plots:  Contour  plots  of  reference  solution  solved  on  extended  domain. 
Full  and  truncated  domains  shown  for  comparison.  Center  Plots:  Contour  plots  of  various 
NRBC  boundary  configurations  using  J  =  1.  Bottom  Plots:  Contour  plots  of  various 
NRBC  boundary  configurations  using  J  =  4 
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Table  5:  Error  as  a  function  of  NRBC  Order  for  open  plane  arbitrary  domain  formu¬ 

lation  using  various  spectral  element  orders  on  square  NRBC  domain.  Gaussian  initial 
condition  and  zero  advection  is  used. 


NRBC  Order 

L  f}  Error 
Linear  Elements 

Error 

Order  2  Elements 

L\  Error 
Order  4  Elements 

J  =  1 

0.13632 

0.13475 

0.13561 

J  =  2 

0.09030 

0.07383 

0.07164 

J  =  3 

0.04483 

0.04403 

0.04293 

J  =  4 

0.03996 

0.03993 

0.03841 

J  =  5 

0.03964 

0.03984 

0.03841 

J  =  10 

0.03948 

0.03957 

0.03841 

J  =  20 

0.03948 

0.03952 

0.03840 

Table  6:  Error  as  a  function  of  NRBC  Order  for  open  plane  arbitrary  domain  formula¬ 

tion  using  various  spectral  element  orders  on  rounded  square  NRBC  domain.  Gaussian 
initial  condition  and  zero  advection  is  used. 


NRBC  Order 

L  jj  Error 
Linear  Elements 

Error 

Order  2  Elements 

L\  Error 
Order  4  Elements 

J  =  1 

0.19805 

0.17845 

0.17331 

J  =  2 

0.08245 

0.07027 

0.06646 

J  =  3 

0.05528 

0.04839 

0.04458 

J  =  4 

0.05246 

0.04599 

0.04358 

J  =  5 

0.05194 

0.04547 

0.04234 

J  =  10 

0.05187 

0.04540 

0.04203 

J  =  20 

0.05186 

0.04540 

0.04201 

C.  ALTERNATIVES 

Thus  far,  using  the  arbitrary  boundary  idea  to  remove  the  problematic  corner  points, 
we  have  not  been  able  to  improve  results  for  the  open  domain  problem.  While  the  arbitrary 
boundary  formulation  does  allow  the  user  to  choose  a  boundary  domain  of  any  shape  (an 
advantage  in  certain  contexts),  the  errors  associated  with  this  formulation  were  on  par  with 
alternatives  presented  in  Chapter  VI.D.  Additionally,  neither  formulation  led  to  exponential 
error  convergence  as  the  order  of  the  NRBC  was  increased.  With  this  in  mind,  we  consider 
an  alternative  boundary  formulation  for  a  circular  domain.  Hagstrom  and  Hariharan  [6] 
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Table  7:  Error  as  a  function  of  NRBC  Order  for  open  plane  arbitrary  domain  formu¬ 

lation  using  various  spectral  element  orders  on  circular  NRBC  domain.  Gaussian  initial 
condition  and  zero  advection  is  used. 


NRBC  Order 

L  jj  Error 
Linear  Elements 

L2)  Error 
Order  2  Elements 

L\  Error 
Order  4  Elements 

J  =  1 

0.26723 

0.26615 

0.25028 

J  =  2 

0.10844 

0.10334 

0.09448 

J  =  3 

0.07584 

0.07211 

0.06451 

J  =  4 

0.07150 

0.06788 

0.06206 

J  =  5 

0.07078 

0.06730 

0.06151 

J  =  10 

0.07070 

0.06720 

0.06142 

J  =  20 

0.07070 

0.06711 

0.06142 

devised  high-order  NRBCs  for  the  standard  time-dependent  two-dimensional  wave  equa¬ 
tion  in  polar  coordinates  implemented  in  a  finite  difference  setting.  This  NRBC  follows 
the  ideas  pioneered  by  Bayliss  and  Turkel  [1]  with  the  exception  that  the  NRBC  condition 
does  not  involve  any  high-order  derivatives  after  introducing  auxiliary  variables.  Huan  and 
Thompson  implemented  the  same  NRBC  in  a  series  of  papers  [54,  55,  56,  57,  58,  59]  in 
a  finite  element  setting.  Here,  we  examine  the  effect  of  this  work  when  examined  with 
high-order  spectral  elements  and  time  integration. 

1.  Hagstrom  Hariharan  Polar  Boundary  Conditions 

The  boundary  condition  devised  by  Hagstrom  and  Hariharan  (hereafter  referred  as 
the  HH  formulation)  provides  a  systematic  approach  for  constructing  boundary  conditions 
for  standard  two-dimensional  wave  equation.  The  condition  is  based  on  the  asymptotic 
series  representation  (which  does  not  converge  at  any  fixed  radius)  for  an  outgoing  solution 
of  the  wave  equation  (in  polar  coordinates) 

1  d2h  d2h  1  dh  1  d2h 

Cq  dt2  dr 2  r  dr  r2  d62  ^ 

Since  the  boundary  condition  is  asymptotic  by  nature,  valid  for  large  radial  distances  -  this 
implies  that  larger  radial  distances  should  provide  better  NRBC  convergence.  Thompson 
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et  al.  make  the  observation  that  . .  for  practical  problems,  truncating  the  asymptotic  ex¬ 
pansion  after  [  J]  terms  provides  solutions  with  errors  well  below  that  of  the  discretization 
error”  [59].  Here,  we  seek  to  significantly  reduce  the  discretization  error  by  employing 
spectral  elements  to  find  the  true  error  convergence  properties  of  the  NRBC.  In  developing 
the  boundary  condition,  Hagstrom  and  Hariharan  construct  a  sequence  of  operators  that 
approximately  annihilate  the  residual  of  the  preceding  element  in  the  sequence,  viewed  as 
a  function  on  the  artificial  boundary.  The  sequence  begins  with  a  first-order  Bayliss-Turkel 
operator  discussed  in  [1],  The  boundary  condition  takes  the  form: 


dh 

dr 

(ftj+i 


0i 


_1_ dh 

c0  dt 


1 

2  r 


h, 


ld4>j  j  {j  ~  \)2  ,  1  d^j-i 

c0  dt  r  '3  4 r2  3  1  4r2  d62 


(VII.7) 
(VII.  8) 


where 


0o  =  2  h  and  0j  =  0. 


At  first  glance,  this  boundary  formulation  suggests  that  we  should  develop  a  “new” 
spectral  element  formulation  for  the  wave  equation  cast  in  polar  coordinates.  If  we  did  this, 
however,  we  would  then  require  a  polar  grid  that  would  introduce  additional  complications 
such  as  the  method  of  dealing  with  the  degenerate  quadrilaterals  that  inevitably  occur  at 
the  center  of  the  grid.  Of  course  there  are  ways  to  overcome  these  obstacles,  but  it  would 
be  much  more  convenient  to  cast  the  problem  in  the  same  framework  already  developed. 
In  other  words,  we  seek  to  implement  this  boundary  condition  (presented  in  polar  form)  in 
our  unstructured  quadrilateral  formulation  of  the  wave  equation  (in  Cartesian  form). 

First,  consider  the  two-dimensional  wave  equation  (same  formulation  as  presented 
in  (II.  29)  with  U  =  V  =  f  =  0) 

S  -  <A2A  =  0.  (VII.9) 
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Multiplying  by  the  test  functions  'I',  and  integrating  over  the  circular  domain  yields  the 


weak  integral  form 


d2h 

t,— r  dn  -  d 


^,;V2h  dn  =  o 


Jn  9t2  Jo. 

Transferring  the  second  order  spatial  derivatives  from  h  to  the  basis  functions  via  integra¬ 
tion  by  parts  and  applying  the  divergence  theorem  to  recast  one  surface  integral  term  as  a 
boundary  integral  gives  us 


T  JO  2 

vl/^  dn  -  c0 


xTin  ■  V h  d,n  +  cl  /  V'kj  •  V h  dn  =  0. 


(VII.  10) 


Of  note  now  is  that  the  boundary  condition  (VII.7)  contains  a  radial  derivatives  of  h  that  on 
the  circle  is  precisely  the  normal  derivative  n  ■  V/i.  This  allows  direct  implementation  of 
the  boundary  condition  into  (VII.  10)  as  follows: 


r  d2h  f  (  1  dh  1  \  f 

Here,  since  on  the  boundary  the  radius  is  fixed,  the  ^  term  may  be  treated  as  a  constant. 

A  similar  weak  form  is  constructed  for  the  boundary  formulation  by  multiplying 
(VII. 8)  by  the  test  functions  Q  and  integrating  over  T  yielding  (after  by  integration  by 
parts): 


I  [  r  d(p3 

co  Jr  ^  9t 


dr  +  3- 

_ 1_  .  end  i  r  QQ  d<j>j- 1 

4r2  *  39  I  start  4 r2  Jr  80  dO 


dT 

dr 


Ci^j+i  dr. 


(VII.  12) 


We  now  use  the  fact  that  the  boundary  is  continuous  and  closed  to  surmise  that  the  endpoint 
evaluation  term  vanishes.  The  formal  problem  statement  is  then:  Find  h  G  V  and  (pj  £  Vp 
where  j  =  1, ...  J  —  1,  such  that  Equations  (VII. 11)  and  (VII. 12)  are  satisfied  V  \Eq  e  V 
and  C;  e  Vr. 
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2. 


Results  for  the  HH  Formulation 


A  series  of  experiments  was  conducted  to  determine  the  effect  of  the  HH  boundary 
condition  for  various  SE  and  NRBC  orders.  Since  the  formulation  is  designed  for  circular 
boundaries,  we  consider  only  circular  boundaries  with  unstructured  grids.  In  each  case,  we 
choose  the  number  of  elements  to  yield  approximately  3,  000  global  points.  Since  the  Gaus¬ 
sian  initial  condition  described  in  (VI. 2)  is  perfectly  symmetric  with  respect  to  the  bound¬ 
ary,  we  introduce  asymmetry  by  adjusting  its  shape  to  yield  a  smooth,  two-dimensional 
“oval-shaped”  initial  condition  with  shape  parameters  ax  =  \,cry  =  |,  further  rotated  by 
an  angle  of  9  —  |.  The  initial  condition  used  here  is: 

h(x,  y,  0)  =  e-{o^+2bxy+a/)  ^  ^  Q)  =  Q  (VII.  13) 


Here,  the  parameters  a,  b,  and  c  are  defined  as  follows: 


cos2  9  sin2  6 

~  2o2  +  2 a2 

x  y 


sin  29  sin  29 
A<72  4cr2 

x  y 


sin2  9  cos2  9 

2  a2  +  2  a2 

x  y 


(VII.  14) 


Again,  the  solution  is  compared  to  one  computed  on  a  larger  domain  allowing  the  wave 
to  propagate  out  of  the  NRBC  domain  but  not  yet  impinge  on  the  non-physical  boundary 
used  to  compute  the  solution  on  the  larger  domain.  Qualitative  results  are  shown  in  Figure 
32  and  quantitative  L\  errors  are  shown  for  various  NRBC  orders  for  SE  orders  up  to  6  in 
Table  8.  No  further  improvement  was  observed  for  SE  orders  above  order  6. 


3.  Adjustments  to  HH  to  Include  Mild  Dispersion 

The  unstructured  grid  representation  of  the  HH  formulation  has  been  demonstrated 
to  significantly  reduce  reflection  caused  by  the  boundary  for  the  standard  wave  equation. 
The  question  now  arises,  can  this  formulation  be  extended  to  include  dispersive  effects  such 
as  Coriolis?  In  [60],  van  Joolen  et  al.  presented  a  method  to  extend  the  HH  formulation 
for  the  standard  wave  equation  under  mild  dispersion.  While  this  formulation  was  well 
grounded  mathematically,  as  far  as  the  author  knows,  it  was  never  implemented.  A  brief 
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(a)  Ref.  Solution  (t  =  1) 


(b)  Ref.  Solution  (t  =  2) 


(c)  Ref.  Solution  ( t  =  3) 


Figure  32:  Open  Domain,  4th  order  spectral  elements  (J  =  4)  using  oblique  Gaussian 
initial  condition  shown  for  t  =  1,  2,  3.  Top  Plots:  Contour  plots  of  reference  solution 
solved  on  extended  domain.  Superimposed  black  circle  indicates  NRBC  domain.  Bottom 
Plots:  Contour  plots  of  various  NRBC  boundary  configurations  using  J  =  4. 

synopsis  of  their  derivation  follows  with  results  presented  for  mild  dispersion  where  f2  = 

0.1. 


We  first  consider  the  KGE  without  advection  (in  polar  coordinates  as  in  the  HH 
derivation): 

1  d2h  r)2h  1  r)h  1  d2h  f2 

-  J—h.  (VII.  15) 


Cq  dt2 


dr  2 


1  dh 
r  dr 


r2  d92  Cn 


As  has  been  previously  discussed,  in  the  geophysical  context,  the  dispersion  parameter  is 
typically  small.  We  assume  here  that 


CqK 


<  1, 


(VII.  16) 
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Table  8:  Tf,  Error  as  a  function  of  NRBC  Order  for  Hagstrom  Hariharan  NRBC  formulation 
using  various  spectral  element  orders  on  the  circular  NRBC  domain.  Oblique  Gaussian 
initial  condition  is  used. 


NRBC  Order 

I2j  Error 
Linear  Elements 

L2}  Error 
Order  2  Elements 

L2,  Error 
Order  4  Elements 

L  jj  Error 
Order  6  Elements 

J  =  1 

0.09310 

0.04772 

0.04555 

0.04485 

J  =  2 

0.03381 

0.00465 

0.00355 

0.00315 

J  =  3 

0.02355 

0.00324 

0.00243 

0.00259 

J  =  4 

0.02217 

0.00305 

0.00236 

0.00201 

J  =  5 

0.02198 

0.00302 

0.00230 

0.00196 

J=  10 

0.02195 

0.00302 

0.00228 

0.00196 

J  =  20 

0.02195 

0.00302 

0.00228 

0.00196 

where  K  is  a  typical  wave  number  appearing  in  the  solution.  Now,  apply  the  Fourier 
transform  to  (VII. 6)  and  (VII.  15)  in  time  to  yield: 


co2 1  d2h 


1  dh 
H — f 


1  d2h 


=  0 


cu2  /2\~  d2h  1  dh  1  d2h  _ 

°  dr2  r  dr  +  r2  dd2 


Wave 

Klein-Gordon 


where  u>  is  the  frequency  and  h  is  the  frequency  domain  representation  of  h.  In  both  cases, 
we  obtain  the  Helmholtz  equation: 

K2h  1  <92/i 

1  dr 2  r  dr  r 2  d6 2 

In  the  non-dispersive  case,  K  —  —  =  K  and  K  =  \  K2  —  ^  in  the  dispersive  case.  In 

CO  y  Cq 

order  to  facilitate  the  conversion  back  to  the  time  domain,  we  now  consider  a  Taylor  series 
approximation  to  the  square  root  term  found  in  the  dispersive  case,  i.e., 

X 

\/l  —  x  —  1  —  -x  +  0(x2). 
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Provided  that  x  is  small,  we  can  truncate  the  0(x 2)  terms.  In  our  case,  from  (VII.  16)  we 
can  reasonably  make  this  assumption  yielding  for  the  dispersive  case: 


K  I  P~  f2  -  f2 

_ —  ,  / 1 - ■- _ ^  1 - - - K  &  K - - - 

K  y  CqA' 2  2c2  IO  2c2  K' 

We  now  see  that  in  the  frequency  domain,  an  equation  valid  in  the  non-dispersive  case  is 
valid  in  the  dispersive  case  if  we  make  the  replacement: 


<viu7) 

We  now  turn  our  attention  to  the  boundary  condition  (VII. 7)  and  (VII. 8)  that  we  Fourier 
transform  in  time  to  yield: 


—iK(j)j  +  -(f) j  — 


-mix  + 


dh  1  * 
dr  2  r  % 


Q-  I)'  1 _ 1  d2<j>j- 1 

4r2  f  ^  4r2  dO2, 


01 

4>j+ 1?  j  =  i)  •  •  •  j  J  i- 


(VII.  18) 
(VII.  19) 


Making  the  substitution  (VII.  17),  we  obtain  the  dispersive  version  of  the  HH  formulation 
in  the  frequency  domain,  i.e., 


.  -  if2  i  dh  1  - 
~lKh  +  MK  +Tr+2rh 


—iK(j)j  + 


if 2 

2c2  K 


V  +  “07  ~ 


(./  )): 

4r2 


A-i 


1  920j_i 
4r2  <902 


01 

0j+ij  7  =  1. 
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Transforming  these  equations  back  into  the  time  domain  results  in  the  final  HH 


boundary  formulation  for  the  KGE: 


i  dh 
c0  dt 


1  d(pj 

c0  dt 


f2  I'1  i 

+  7TT  /  0i(r)dr+-0i 
zco  Jo  r 

" - V - ' 

n(t) 


f 


dh 


—  /  h(T)dT+—  +  —  h 


2  Ci 


o  Jo 


U-i) 

4  r2 


m(t) 

2 


1 


<9r  2r 


1 

4r2  (96*2 


0i 

0j+l 


(VII.  20) 


(VII.21) 


where 


j  =  1,  •  •  • ,  J  ~  1,  0o  =  2 h  and  0j  =  0. 


It  should  be  noted  that  van  Joolen  et  al.  [60]  show  how  m(f)  and  n(f)  can  be  calculated 
in  each  time- step  to  keep  the  boundary  condition  local  in  time  without  having  to  store  and 
operate  on  the  history  of  the  solution.  For  this  analysis,  a  simple  trapezoidal  approximation 
was  used  to  approximate  the  integral. 

The  weak  form  of  the  formulation  is  now  constructed.  We  consider  the  KGE  in  its 
general  form: 

-  CqV2/i  +  f-h  =  0. 

Multiplying  by  the  test  functions  T,  and  integrating  over  the  circular  domain  yields  the 
weak  integral  form 

/  (in  -  c2  [  't>l\72hdtt  +  f2  I  Vihdn  =  0. 

Jn  dt  Jq  Jq 


Transferring  the  second  order  spatial  derivatives  from  h  to  the  basis  functions  via  integra¬ 
tion  by  parts  and  applying  the  divergence  theorem  to  recast  one  surface  integral  term  as  a 
boundary  integral  gives  us 


d2h 

*‘Wdn~c2° 


■  V/i  dn  +  cl  /  VTq  •  Vh  dfl  +  f2  /  dfl  —  0.  (VII. 22) 
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Of  note  now  is  that  the  boundary  condition  (VII.20)  contains  a  radial  derivatives  of  h  that 
on  the  circle  is  precisely  the  normal  derivative  ft  ■  Vh.  This  allows  direct  implementation 
of  the  boundary  condition  into  (VII.22)  as  follows: 


+c0 


1  dh  1  f 2  \ 

W*  ■  Vh  dn  +  f  I  ^ih  dQ  =  0. 
i  Jn 


(VII.23) 


Here,  since  on  the  boundary  the  radius  is  fixed,  the  ^  term  may  be  treated  as  a  constant. 

A  similar  weak  form  is  constructed  for  the  boundary  formulation  by  multiplying 
(VII.21)  by  the  test  functions  Q  and  integrating  over  T  yielding  (after  by  integration  by 
parts): 


1 

c0 


<90 


/ 


Ci-^dT  +  f-  /  Qn(t)dT+J-  /  Ci^dT 


(J-IY 


dt 


2  cL  . , 
4r2 


<90j_i  iend 


dd 


r  _  ^ 

1 

start  4  r2 


3  4r2 

dQ  d(t>j- 1 


dd  dd 


dT  = 


QQj—i  dr 

r 

Q,0j  \- 1  dr. 


(VII.  24) 


We  now  use  the  fact  that  the  boundary  is  continuous  and  closed  to  surmise  that  the  endpoint 
evaluation  term  vanishes.  The  formal  problem  statement  is  then:  Find  h  G  V  and  <p:]  e  Vr 
where  j  =  1, ...  J  —  1,  such  that  Equations  (VII.23)  and  (VII. 24)  are  satisfied  V  Tq  e  V 
and  Q  e  Vr. 


4.  Results  for  HH  with  Dispersion 

A  series  of  experiments  was  conducted  to  determine  the  effect  of  the  HH  boundary 
condition  extended  to  include  mild  dispersion  for  various  SE  and  NRBC  orders.  The  set-up 
is  identical  to  the  experiments  without  dispersion,  except  the  dispersion  parameter  is  set  to 
f2  =  0.1.  Qualitative  results  are  shown  in  Figure  33  and  quantitative  Ip  errors  are  shown 
for  various  NRBC  orders  for  SE  orders  up  to  6  in  Table  9.  As  in  the  non-dispersive  case, 
no  improvement  was  observed  for  SE  orders  above  order  6. 
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(a)  Ref.  Solution  (t  =  1) 


(b)  Ref.  Solution  (t  =  2) 


(c)  Ref.  Solution  ( t  =  3) 


(d)  NRBC  J  =  4  (t  =  1) 


-0.1  -0.05  0  0.05  0.1 


(e)  NRBC  J  =  4  (t  =  2) 


(f)  NRBC  J  =  4  (f  =  3) 


Figure  33:  Open  Domain,  Ath  order  spectral  elements  (J  =  4)  using  oblique  Gaussian 
initial  condition  shown  for  t  =  1,  2,  3  under  dispersion  f2  =  0.1.  Top  Plots:  Contour  plots 
of  reference  solution  solved  on  extended  domain.  Superimposed  black  circle  indicates 
NRBC  domain.  Bottom  Plots:  Contour  plots  of  various  NRBC  boundary  configurations 
using  J  =  4. 
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Table  9:  Error  as  a  function  of  NRBC  Order  for  Hagstrom  Hariharan  NRBC  formulation 

using  various  spectral  element  orders  on  the  circular  NRBC  domain.  Oblique  Gaussian 
initial  condition  is  used  with  dispersion  parameter  set  to  / 2  =  0.1. 


NRBC  Order 

L  f,  Error 
Linear  Elements 

L\  Error 
Order  2  Elements 

L’o  Error 
Order  4  Elements 

7.  if  j  Error 
Order  6  Elements 

J  =  1 

0.07290 

0.03555 

0.03369 

0.03293 

J  =  2 

0.02684 

0.00371 

0.00283 

0.00248 

J  =  3 

0.01869 

0.00258 

0.00192 

0.00204 

J  =  4 

0.01759 

0.00243 

0.00186 

0.00157 

J  =  5 

0.01744 

0.00240 

0.00181 

0.00154 

J=  10 

0.01742 

0.00240 

0.00180 

0.00153 

J  =  20 

0.01742 

0.00240 

0.00180 

0.00153 
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VIII.  CONCLUSIONS  AND  AREAS  FOR  FUTURE  RESEARCH 


In  this  dissertation,  we  considered  a  reduced  form  of  the  shallow  water  equations  in 
various  (semi-infinite  and  infinite  channels  as  well  as  open  domain)  configurations.  Using 
the  Givoli-Neta  auxiliary  variable  formulation  of  the  Higdon  non-reflecting  boundary  con¬ 
ditions,  we  truncated  the  original  infinite  domain  and  developed  the  boundary  conditions 
specific  to  the  problem  at  hand.  Using  a  high-order  approach  to  the  spatial  discretization 
(spectral  elements),  time  integration  (high-order  Runge-Kutta)  in  concert  with  high-order 
boundary  treatment,  we  showed  exponential  convergence  to  the  reference  solution  in  chan¬ 
nel  configurations.  These  results  suggest  a  balanced  approach  to  dealing  with  truncation 
errors  -  namely,  to  make  improvements  in  all  components  of  the  problem  to  see  improved 
accuracy. 

In  open  domain  problems,  we  considered  various  ways  to  handle  corner  compatibil¬ 
ity  concerns  when  using  the  Givoli-Neta  auxiliary  variable  formulation.  Using  a  physical 
argument  that  the  auxiliary  variables  themselves  should  be  non-reflecting  at  a  boundary, 
we  formulated  a  spectral  element  formulation  that  yielded  stable  solutions  (even  for  long 
term  time  integrations)  using  first  order  NRBCs  for  the  auxiliary  variables.  This  formula¬ 
tion  showed  significant  improvement  from  the  first  order  ( J  —  1  Sommerfeld  condition) 
to  higher  order  J,  although  the  improvement  was  far  from  exponential.  Besides  the  low 
order  method  of  handling  the  auxiliary  variable  boundaries,  the  “node-splitting”  method  of 
handling  double-counting  corner  nodes  turned  out  to  be  convergence  limiting. 

Recognizing  that  any  formulation  that  included  corners  would  be  problematic,  we 
sought  a  boundary  condition  that  would  be  valid  for  an  arbitrarily  shaped  domain.  The 
Givoli-Neta  formulation  was  shown  to  have  insurmountable  implementation  issues  without 
the  simplifying  assumption  of  small  curvature  on  the  boundary.  When  the  small  curvature 
assumption  was  made,  the  formulation  was  shown  to  have  stable,  improved  results  from 
the  first  order  ( J  =  1)  Sommerfeld  condition.  It  was  clear,  however,  that  there  are  trade¬ 
offs  associated  between  accurately  representing  the  Givoli-Neta  formulation  and  removing 
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the  problematic  corners.  Specifically,  using  a  square  domain  (where  the  small  curvature 
assumption  affects  only  the  four  comer  points)  errors  were  shown  to  be  improved  over 
alternative  domains  (rounded  square  and  circle). 

The  final  experiments  conducted  considered  a  boundary  condition  originally  de¬ 
vised  by  Hagstrom  and  Hariharan  and  extended  to  the  dispersive  wave  equation  by  van 
Joolen  et  al.  This  boundary  condition  is  based  on  the  asymptotic  series  representation 
of  the  wave  equation  in  polar  coordinates,  valid  for  large  radial  distances.  Results  were 
improved  over  the  alternative  arbitrary  domain  formulations  although  valid  only  for  large 
radial  distances  and  restricted  to  circular  boundary  domains  in  this  analysis. 

This  research  has  demonstrated  exponential  convergence  in  channel  experiments, 
only  hypothesized  in  previous  low-order  settings.  Additionally,  it  has  developed  several 
alternatives  to  handling  open  domains  which  improve  performance  to  first-order  NRBC 
schemes  at  a  very  moderate  computational  cost.  What  remains  is  to  extend  this  high-order 
numerical  formulation  to  more  complex  linear  and  non-linear  systems  of  fluid  motion  such 
as  the  Euler  equations.  Additionally,  better  alternatives  to  dealing  with  corner  compatibility 
concerns  remains  an  open  problem. 
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APPENDIX  A.  DEPTH  INTEGRATING  THE  CONTINUITY 

EQUATION 


In  order  to  arrive  at  the  final  form  of  the  shallow  water  equations,  we  depth  inte¬ 
grated  our  shallow  water  continuity  equation.  The  details  follow.  Given 


0  =  V  ■  u 


we  integrate  in  z 


0=1  V  ■  u dz 

J-hB 

fh  ( du  dv.  1 

=  .L  U  +  aS,l<fa  +  wU-w|-^ 


(A.l) 


Since  both  h  and  h  B  depend  on  x  and  y  (and  t  for  h ),  we  apply  the  Leibniz  integral  rule, 
which  allows  us  to  write: 


d_ 

dx 

d 


r-z=h 

r-z=h 

du  , 

dh 

/  udz  = 

— —  dz  +  u 

— - u 

'  z=—hB 

J  z=-hB 

dx 

hdx 

pz=h 

f‘Z=h 

dv 

dh 

/  vdz  = 

-z—dz  +  v 

— - V 

J z=—hB 

J  z=-hB 

dx 

hdx 

d{-hB) 

hB  dx 

d(-hB) 

hB  dx 


Substituting  these  into  (A.l)  we  have 


d  fh  , 

0  =  —  uaz  — 
dx  J — foB 

9  fh 

+  —  /  vdz  - 

dy  J~hB 


u 


dh 

z=h  dx 

dh 

‘=hdy 


—  u 


d(-hB 

t=—hB  dx 

d  (— hB ) 

z=-hB  dy 


w 


w 


z=h 


z=—hf 
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Now,  using  the  appropriate  boundary  conditions 


w(x,y,-hB) 
w(x,  y,  h,  t) 


d  hB  dhB 

-u— - V- 


dx 


dy 


dh  dh  dh 
si  +  udx+v8y 


we  substitute  to  find 


d  fh  ,  d  fh  , 

0  =  —  /  uaz  +  —  /  vdz 
9x  J-hB  dy  J-hB 

dh 

+  di 


-u 


+u 


dh 

z=h  dx 

dh 

z=h  dx 


dh 
*=h  dy 
dh 
t=hdy 


u 


u 


dhB 

z=—hB  dx 

dhB 

z=-hB  dx 


+  V 


dhB 

z—  hB  dy 

dhB 

z=-hB  dy 


Which  simplifies  to 


dh  d  fh  d  fh 

—  +  —  /  udz  +  —  /  vdz. 
dt  dx  J_hB  dy  J^hB 


The  construction  of  the  shallow  water  model  as  shown  in  Figure  3  has  H  —  h  +  hB  where 
H  is  the  depth  of  the  fluid.  Since  a  previous  argument  showed  that  u  and  v  are  independent 
of  depth,  we  are  left  with  our  final,  depth  integrated  continuity  equation 


dh 

dt 


d_ 

dx 


(Hu)  + 


d_ 

dy 


(. Hv )  =  0 
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APPENDIX  B.  LINEARIZING  THE  SHALLOW  WATER 

EQUATIONS 


The  non-linear  version  of  the  shallow  water  equations  are: 

dtu  +  udxu  +  vdyu  —  fv  =  —g  dxh 
dtv  +  udxv  +  vdyV  +  fu  —  —g  dyh 
dth  +  dx  (Hu)  +  dy  (. Hv )  =  0. 

We  wish  to  find  a  linear  version  of  these  equations.  Suppose  the  bottom  topography  is  flat 
such  that  hB  is  constant  and  u  and  v  can  be  described  by  a  constant  mean  term  and  a  small 
0(5)  deviation  from  that  value,  i.e., 

u  =  U  +  u*  v  —  V  +  v*  H  =  Jib  +  h 

To  be  clear,  U  and  V  are  the  mean  velocities  and  hB  is  the  mean  water  depth.  We  now 
make  these  perturbation  substitutions. 

d t(U  +  u*)  +  (U  +  u*)dx(U  +  u*)  +  (V  +  v*)dy(U  +  u*)  -  f(V  +  v*) 

dt(V  +  v*)  +  (U  +  u*)dx(V  +  v*)  +  (V  +  v*)dy{V  +  v*)  +  f(U  +  u*) 

d th  +  dx  (( hB  +  h)  (U  +  u*))  +  dy  (( hB  +  h)  (V  +  v*)) 

Now,  recalling  that  U,  V  and  hB  are  constants,  we  simplify  to  find: 

dtu*  +  (U  +  u*)dxu*  +  (V  +  v*)dyu*  —  f(V  +  v*)  =  —g  dxh 

3tv*  +  (U  +  u*)dxv*  +  (V  +  v*)dyv*  +  f(U  +  u*)  =  -g  dyh 

dth  +  hB  (dxu*  +  dyv*)  +  dxh  (U  +  u*)  +  dyh  (V  +  v*)  +  h  (dxu*  +  dyv *)  =  0. 


=  ~g  dxh 

=  -g  dyh 
=  0. 
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If  we  expand  each  term,  the  result  is: 


d tu*  +  Udxu*  +  u*dxu*  +  VdyU*  +  v*dyu*  —  f(V  +  v*)  =  —g  dxh 
dtv*  +  Udxv*  +  u*  dxv*  +  Vdyv*  +  v*dyv*  +  f{U  +  u*)  =  -g  dyh 
dth  +  hsdxU*  +  hsdyV*  +  Udxh  +  u*dxh  +  Vdyh  +  v*dyh  +  hdxu*  +  hdyv*  =  0. 

We  now  neglect  any  terms  of  0(d2)  to  arrive  at  our  final  form  of  the  linearized  shallow 
water  equations: 


dtu*  +  Udxu*  +  Vdyu*  —  f(V  +  v*)  =  —g  dxh 
dtv*  +  Udxv*  +  Vdyv*  +  f{U  +  u*)  =  -g  dyh 
dth  +  U dxh  +  Vdyh  +  hs  ( dxu *  +  dyv *)  =  0. 
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APPENDIX  C.  ADJUSTING  HIGDON’S  CONDITION  FOR 

ADVECTION 


Suppose  that  we  have  a  wave  that  moves  according  to 

(dt  +  Udx)2  h  -  c20dlh  =  0.  (C.l) 

This  equation  can  be  “factored”  as  follows: 

0  =  (^  dt  +  (U  —  Co)  dx^j  (dt  +  (U  +  Co)  dx^j  h 

Following  a  standard  method  of  characteristics  derivation,  we  define  functions  w  and  v  as 

w(x,  t )  =  d th  +  (U  —  c0)dxh 
v(x,  t )  =  dth  +  (U  +  c0)dxh. 

It  can  be  verified  that  both 

d tw  +  (U  —  c0)  dxw  =  0 
dtv  +  (U  +  c0)  dxv  =  0 

satisfy  (C.l)  exactly.  If  we  then  solve  w  along  its  characteristic  x  =  (U  —  c0)  t  +  x0  and 

v  along  its  characteristic  x  =  (U  +  c0)  t  +  x0,  with  initial  data  w(x0,  0)  =  P(x0)  and 

v(x0,  0)  =  Q(x o),  then  the  solutions  for  w  and  v  are  respectively 

w(x,t )  =  P^x  —  (U  —  c0)t ^  =  dth  +  (U  —  c0)dxh  (C.2) 

v(x,  t)  =  Q^x  -  (U  +  c0)tj  =  dth  +  (U  +  c0)dxh  (C.3) 
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We  subtract  (C.3)  from  (C.2)  to  find 

p[x  —  (U  —  c0)t  )  =  d th  +  (U  +  c0)dxh 
-  q(x-  (U  +  c0)t  )  =  d th  +  (U  -  c0)dxh 

P(x-(U-  c0)t )  -  Q  [x  -  (U  +  c0)t)  =  2codxh.  (C.4) 

Further,  if  we  combine  (C.2)  and  (C.3)  as 

(U  -  co)  P(x  -{U-  c0)t)  =(U-  c0)  (dth  +  (U  +  c0)dxh ) 
—  (U  +  Co)  q[x  —  (U  +  co)f)  =  (U  +  Co)  [dth  +  (U  —  co)dxh^j 

(U  -  c0)  p[x-(U-  c0)t)  ~{U  +  c0)  q(x-(U  +  c0)tj  =  -2c0dth.  (C.5) 

This  implies  that  the  solution  takes  the  form  h(x,  t)  —  F  (U +c0)t  )  +G  (U— c0)t  ) . 

Here,  F  and  G  are  arbitrary  functions  of  the  initial  data.  To  see  this,  consider 

dth  =  -{U  +  co)F'{x  -  ( U  +  c0)t)  -  (U-  c0)G'(x  -  (U  -  c0)t)  (C.6) 

dxh  =  F'[x-  (U  +  c0)tj  +G'(x-(U-  c0)t) .  (C.7) 

Equating  coefficients  with  (C.4)  and  (C.5)  this  yields  a  relation  between  the  initial  data  and 
the  functions  F  and  G 

p[x  —■  {U  —  co)*)  =  2c0G"  (x-(U-  co)*)  (C.8) 

q(x-(U  +  co)*)  =  -2c0F' (x  -  (U  +  co)*) •  (C.9) 

This  solution  can  be  rewritten  as 

/i(x, * )  =  F (x  -  (c o  +  U)t )  +  G (x  +  (c0  -  U)tj  (C.  10) 
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with  the  interpretation  that  the  general  solution  is  the  sum  of  F  yx  —  (co  +  U)tJ ,  a  wave  of 
fixed  shape  moving  to  the  right  with  velocity  (c0  +  U)  and  G  +  (c0  —  U)tj  a  wave  of 
fixed  shape  moving  to  the  left  with  velocity  (c0  —  U). 
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APPENDIX  D.  NORMAL  TO  TANGENTIAL  DERIVATIVE 
TRANSFORMATION  FOR  EASTERN  NRBC 


The  function  h  satisfies  the  dispersive,  advective  wave  equation  (11.29)  in  D.  Since 
the  function  di  is  obtained  by  applying  the  linear  (constant  coefficient)  operator  [dx  +  ^rdt 
to  h,  it  is  can  be  shown  that  <fii  should  also  satisfy  the  same  equation  in  D9 .  Further,  since 
(f)j  is  obtained  by  applying  the  same  linear  operator  j  —  1  times  to  0\,  the  functions  <P:] 
should  satisfy  an  equation  like  (11.29),  namely, 

(dtt  +  (  U 2  —  Cg)  dXX  +  (V"  —  Cq)  dyy  + 

2  Udxt  +  2  Vdyt  +  2  UVdxy  +  /2)  <j)j  =  0  (D.l) 

Now,  use  the  following  identities: 

dxx4>j  —  (dx  —  T7  j  (dx  +  —  dt  j  4>j  + 

V  L'j+ 1  J  \  W+i  J  L’j+ i 

^xi0j  dt  ( dx (pj) 
dxy0j  dy  ( dx  0-j ) 

and 

dx  + -^dt^  (pj-i  =  (pj  j  =  l,...,J.  (D.5) 

9Here  we  must  use  the  assumption  that  Co  and  /  are  constants.  By  applying  the  differential  operator  to 
(11.29),  computing  each  of  the  (pj  derivatives  present  in  (III.  17)  using  the  differential  operator  and  simplifying, 
a  simple  induction  argument  shows  that  the  (pj’s  must  satisfy  (III.  17) 


(D.2) 

(D.3) 

(D.4) 
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Now,  if  we  substitute  (D.2)  -  (D.4)  into  (D.l),  and  replace  j  with  j  —  1  everywhere 
yields,  for  j  =  1, . . . ,  J 

fa- 1  +  (U2  ~  Co)  (dx  ~  c'dt'j  fa- 1  +  Q2fa-i  +  ( Y 2  -  co)  + 

2f/  dt  (dx(f)j- 1)  +  2E  dyt4>j- 1  +  2 UVdy  (dx<f>j- 1)  +  f24>j- 1  =  0  (D.6) 

From  this  and  (D.5)  one  gets,  for  j  =  1, . . . ,  J 

fa-i  +  (^2  _  co)  $3  +  (vfa-i  +  (^2  _  co)  fayfa- 1+ 

2 Udt  (dx(j>j—i)  +  2Vdyt4>j~i  +  2UVdy  (dx<f)j-\)  +  f2(f>j- 1  =  0  (D.7) 

Now,  we  shift  indices  on  (D.5)  and  multiply  both  sides  by  (f/2  —  Cq)  as: 

(172  -  c2)  ^  4>j  =  ( U 2  -  c20)  03+l  j  =  0, . . . ,  J  -  1.  (D.8) 

Subtract  (D.7)  from  (D.8) 

~fa-i  +  ( u 2  -  Cq)  fa  -  (^2  _  co)  £2^—1  _  -  co)  fayfa- 1- 

2f/9t  (ax0J-_1)  -  2Vdvtfa.1  -  2 UVdy  (dxfa^)  -  ffa. 1  =  ((72  -  c2)  0J+1  (D.9) 

Now,  consider  (D.5).  Expanding  and  solving  for  dx<pj_ i,  we  get: 

dxfa-i  =  fa  -  ^rfa-i  j  —  1, . . . ,  J  .  (D.10) 
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Substitute  (D.10)  into  (D.9) 


~h- 1  +  (^2  ~  co)  ^  ~  (^2  _  co) 

( v 2  -  Co)  dyyfa- 1  -  2Udt  (^j)j  -  -  2Vdyt(j)j-i-  ( D .11) 

2UVdy  (t,  -  i  =  (t/2  -  c2)  0J+1 

Simplify: 


U2  —  co 


2C/1/ 


2K  ^■_1  -  (V2  -  c2)  0".,+ 


(V  -  Co)  (2-  +  p)-)  -  2 uj  4,  -  2UV4>'j  -  /Vj-,  =  (l r2  -  cl)  <fe+i  (D.12) 


In  (D.12)  and  elsewhere,  a  prime  indicates  differentiation  with  respect  to  y  along  T e,  i.e. 
the  tangential  derivative  along  those  boundaries. 
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APPENDIX  E.  METRIC  TERMS  DERIVATION 


To  facilitate  interpolation  and  integration  required  for  the  SE  method,  we  transform 
all  terms  of  the  weak  integral  form  in  physical  space  x  =  (x.  y)1  to  a  canonical  space  £  = 
(£,  r])T.  This  nonsingular  mapping  assumes  x  =  x  (£,  rj)  and  conversely  ^  =  $,  (x,  y). 

A.  DERIVATION  OF  METRIC  TERMS 

Using  the  chain  rule,  we  find 

,  dx  dx 
dx  =  —d£  +  —dr/, 

which  can  be  written  in  matrix  form 


Here  Je  is  the  transformation  Jacobian  with  associated  determinant  \Je\  defined  as 


dx 

dx 

d(, 

dy 

dy 

dy 

at. 

dy 

dx  dy  dy  dx 
d£  drj  d£  drj ' 


Similarly,  we  can  find  the  inverse  transformation 


d£  =  it dx  +  ^ dy 
dx  dy 


and  write  the  derivatives  of  £(x,  y)  in  matrix  form  as 
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where 


Te 


d£  d£ 

dx  dy 
dy  dy 
dx  dy 


(E.4) 


is  the  Jacobian  of  the  inverse  transformation.  Since  the  Jacobian  described  by  (E.2)  is 
non-singular,  we  can  write  the  transformation  described  by  (E.l)  as 


drj 


Here  ( Je)  1  is  the  standard  matrix  inverse  defined  by 


(JT1  = 


dy  _ dx ■ 

dy  dy 

dy  dx 


(E.5) 


(E.6) 


We  now  note  that  the  formulations  described  by  (E.3)  and  (E.5)  are  identical,  therefore, 
equating  coefficients  from  the  two  Jacobian  terms  (E.4)  and  (E.6)  yields  the  metric  terms 


<9£  1  dy  <9£  1  dx  drj  1  dy  dr)  1  dx 

dx  \Je\dr)'  dy  \Je\dr)'  dx  \Je\  <9£’  dy  \Je\  dC, 

This  formulation  is  convenient  since  all  metric  terms  are  defined  in  terms  of  terms  that  are 
readily  calculated  via  basis  function  expansions  of  cc(£,  rj). 


B.  CONSEQUENCES  OF  QUADRILATERAL  GRID  DEGENERATION 


Given  the  discussion  of  metric  terms,  we  see  that  there  is  a  common  term  in  that  has 
the  potential  to  cause  numerical  instability  -  namely  hL.  If  the  elemental  Jacobian  tends  to 
zero,  then  all  metric  terms  associated  with  that  Jacobian  will  tend  toward  infinity.  In  fact, 
if  the  element  Jacobian  is  small  at  certain  points  on  the  global  domain  compared  with  other 
locations  on  the  global  domain,  experiments  in  this  dissertation  have  shown  that  they  tend 
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to  corrupt  the  entire  solution.  The  question  then  arises,  what  element  geometries  have  the 
potential  to  cause  numerical  instabilities? 

Since  the  global  problem  is  discretized  into  smaller  elements,  this  question  must  be 
answered  in  the  context  of  grid  generation.  There  is  general  agreement  in  literature  [45, 
61,  62,  63,  64,  65,  66]  concerning  quadrilateral  grid  generation  that  convex  elements  with 
maximum  internal  angles  ~  135°  constitute  a  quality  mesh.  Li  et  al.  in  [67]  describe  pro¬ 
cedures  to  adjust  quadrilateral  basis  functions  to  deal  with  elements  that  have  large  interior 
angles  or,  in  fact  completely  degenerate  into  triangles,  although  with  these  adjustments, 
computational  overhead  is  increased  to  deal  with  alternate  canonical  geometries.  In  this 
analysis,  we  seek  a  quality  all  quadrilateral  mesh. 

To  demonstrate  how  a  poorly  generated  mesh  can  taint  a  solution,  we  consider  an 
extreme  example  where  a  quadrilateral  element  is  degenerated  into  a  triangle  as  shown  in 
Figure  34.  In  this  case,  one  of  the  internal  angles  of  the  “ quadrilateral ”  is  180°. 


Figure  34:  Degenerate  quadrilateral  mapped  to  a  canonical  reference  element. 
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Now,  consider  the  linear  basis  function  expansion  (similar  argument  for  higher  or¬ 
der  expansions)  of  the  physical  coordinates  and  their  derivatives 

4 

=  ^2^k(t,  v)  xk 

k=  1 

)  =  y'  ) 

where  the  linear  Lagrange  basis  functions  are  defined  as 

V’i  (f ,  v)  =  I  (1  -  0  (!  -  v) ,  ^2  (f ,  v)  =  \  (1  +  0  (!  -  v)  r 

^3  (f ,  v)  =  \  (!  -  0  (!  +  v) ,  fa  (f ,  v)  =  \  (1  +  0  (!  +  v)  ■ 

We  consider  the  degenerate  vertex  located  at  (£,  rf)  =  (1, 1)  and  compute  each  term 
required  for  computation  of  the  Jacobian. 

^  (1, 1)  =  2  (-*3  +  XA) ,  ^  (1, 1)  =  -  (~x2  +  x4) , 

dy  1  dv  1 

^  (1,  1)  =  2  (-1/3  +  1/4) ,  ^  (1, 1)  =  2  ("Ifc  +  1/4)  • 

Using  the  fact  that  (x2,  y2) ,  (x3,  y3) ,  (x4,  y4)  are  collinear,  we  put  x4  in  terms  of  the 
other  two  points  yielding 

*4=  (y4-y3  +  x3(^^-))  f^V  (E.8) 

V  \x3-x2JJ  \y3-y2J 

Computing  the  determinant  of  the  Jacobian  (E.2)  and  simplifying  using  (E.8)  yields 

\Je\  =  {~x3  +  xA)  (— y2  +  Vi)  ~  (-J/3  +  Va)  {~x2  +  x4) )  =  0  (E.9) 

This  degenerate  quadrilateral  will  have  infinite  metric  terms.  Even  if  the  quadrilateral  el¬ 

ement  does  not  completely  degenerate  into  a  triangle,  but  has  a  large  angle  -  the  metric 
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terms  will  be  very  large  in  comparison  to  other  element  metric  terms  and  will  have  a  desta¬ 
bilizing  effect.  Therefore,  all  meshes  in  this  analysis  were  generated  with  internal  angles 
less  than  135°  choosing  to  add  additional  elements  (and  degrees  of  freedom)  to  ensure  that 
this  happens. 
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APPENDIX  F.  N ON -DIMENSIONALIZ ATION  OF  THE  KGE 


Consider  the  KGE 


(dt  +  Udx  +  Vdyf  h  -  CqV2/t  +  fh  =  0 


that  we  wish  to  study  in  a  non-dimensional  context.  For  simplicity  in  derivation,  we  assume 
that  there  is  no  advection  (U  —  V  —  0),  yielding: 

^-4v*h+n= o. 


Now,  as  outlined  in  [51],  we  examine  typical  scales  of  motion  in  the  ocean  so  as  to  recast 
the  problem  in  a  dimensionless  way  (can  substitute  typical  scales  for  atmosphere  or  other 
medium  as  well  with  the  same  process  that  follows).  For  this  analysis,  the  length  scales 
were  chosen  0(100  km),  vertical  depth  scales  0(100  m),  scales  for  h  0(1  m)  and  the 
dispersion  parameter  /  for  Coriolis  O  (10~4s_1).  Given  these  choices,  we  know  from  the 
discussion  in  Chapter  II  that 

/  m  \  m2 

cn  =  9^b  =  ( 10—  )  (100m)  =  1000  — 

V  s2/  sz 


that  makes  our  specific  problem: 


d2h 
dt 2 


V2h  + 


0. 


Now,  we  follow  the  details  as  outlined  in  [68]  to  scale  out  any  dimensions.  In  particular, 
we  define  the  following: 

x  m  2/m  —  h  m 

x  =  — ^  —  y  =  — -  —  h  — - 

105  m  105  m  1  m 
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where  the  typical  length  scale  is  100  km  (105  m).  Now,  we  note  via  the  chain  rule  that: 


dh  dh  dx 

dx  dx  dx 

d2h  d  (dh \ 

dh  1 

dx  105m 
!  9 

(dh\ 

,  1  d  / 

dx2  dx  ' 

l  dx  j 

105m  dx  ' 

{dx  ) 

105m  dx  \ 

d2li 


1010m2  dx2 


similarly, 


d2h 


d2h 


dy 2  1010m2  dy2 


Using  this  information  in  (F)  yields: 

d2h  1000m2 


d  t2 
d2li 


1010m2 


(v2ft) 


10~  I  „ 

H - tt-  h  —  0 


1  v2h  +  ^h  =  o 


dt2  107s2 

>107  s2^  -  \72h  +  0.lh  =  0 
dt 2 


Now,  we  remove  the  dimensions  from  our  variable  h  to  get 


107  s2^  -  V2h  +  0.1h  =  0. 
dt 2 


Letting  t  =  |(|[-  ^  and  noting  via  a  similar  argument  as  above  that 

d2h  _  1  d2h 

dt 2  107  s2  dt2 

we  arrive  at  our  final,  non-dimensional  form  of  the  Klein-Gordon  Equation, 

d2h 


dt 2 


V2/r  +  0.1  h  =  0 


where  t  =  (103  -5  s)  t,  h  =  (1  m)  h,  x  =  (10°  m)  x  and  y  =  (105  m)  y. 
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APPENDIX  G.  AUXILIARY  VARIABLE  LORMULATIONS  LOR 
WESTERN,  NORTHERN  AND  SOUTHERN  BOUNDARIES 


For  configurations  studied  in  this  dissertation,  G-N  auxiliary  variable  formulations 
are  required  for  boundaries  other  than  T e,  explicitly  derived  in  Chapter  III.  What  follows 
here  are  the  details  of  the  formulation  for  each  of  the  other  boundaries. 

A.  FORMULATION  FOR  THE  WESTERN  BOUNDARY 


We  begin  by  stating  the  Higdon  boundary  condition  for  Tw  given  by: 


Hj : 


0  on  r  iy  . 


(G.l) 


When  imposed  as  a  boundary  condition  on  Tw,  we  can  recast  this  formulation  in  terms  of 
auxiliary  variables  as  outlined  in  Chapter  III  equivalent  to  the  single  boundary  condition 
(G.l)  as: 


where:  o0  =  h 


j  =  (G.2) 


4>j  =  o 


This  set  of  conditions  involves  only  first-order  derivatives.  However,  due  to  the  appear¬ 
ance  of  the  ^-derivative  in  (G.2),  one  cannot  discretize  the  o3  on  the  boundary  Tw  alone. 
Therefore,  we  shall  manipulate  (G.2)  in  order  to  get  rid  of  the  rc-derivative. 

The  function  h  satisfies  the  dispersive,  advective  wave  equation  (11.29)  in  I).  Since 
the  function  0\  is  obtained  by  applying  the  linear  (constant  coefficient)  operator  (dx  -  ^rdt 
to  h,  it  is  clear  that  <i>\  should  also  satisfy  the  same  equation  in  1).  Further,  since  03  is  ob¬ 
tained  by  applying  the  same  linear  operator  j  —  1  times  to  ©i .  the  functions  4>j  should  satisfy 
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an  equation  like  (11.29),  namely, 


(dtt  +  (U~  —  Cq)  dxx  +  —  Cq)  dyy+ 

2  Udxt  +  2U  dyt  +  2I/UV,  +  /2)  y  =  0  (G.3) 

Using  the  following  identities: 

dxxrf’j  —  (dx  +  -q  dt  \  (dx  —  -q  dtj  4>j  +  0j 

\  W+i  /  V  V'+i  /  Ly+i 

()x  l0j  ((/  ( 0j  ) 

^xy0j  (  A/  (  0r  0j  ) 


and  combining  with  (G.2)  allows  us  to  write  (G.3)  as: 


/  2C7  U2-c2\ 

V  Q  q  ) 


h-1  -  [y  +  2v)  -  (l-22  -  c2) 

(4  +  V)  +  2C/)  -  2UV^  -  /Vi- 1  = 


((72  -  c2)  0,+1 


for  j  =  l,...,J-l  (G.4) 


In  (G.4)  and  elsewhere,  a  prime  indicates  differentiation  with  respect  to  y  along  IV,  be., 
the  tangential  derivative  along  Tw.  As  desired,  the  new  boundary  condition  (G.4)  does  not 
involve  ^-derivatives.  In  addition,  there  are  no  high -y  or  t  derivatives  beyond  second  order. 
The  new  formulation  of  the  Jth- order  NRBC  on  Yw  can  be  summarized  as  follows: 

Poh  +  dxh  =  , 

aj4>j- 1  +  1  —  +  /3j4>j  —  70)  —  /20j— 1  —  Ax0j+i  (G.5) 

0o  =  h  <t>j  =  0 
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where 


A)  — 


aj 


2  U  U2  -  eg 

c,  q 


Pi 


-2U, 


Kj  = 


2UV 

~C~ 


—  2V, 


7  =  2UV, 


\y  =  V2~  Cg, 

\x  =  u2  -  cl 


B.  FORMULATION  FOR  THE  NORTHERN  BOUNDARY 


The  Higdon  boundary  condition  for  T N  is  given  by: 


Hj  : 


(G.6) 


When  imposed  as  a  boundary  condition  on  T n,  we  can  recast  this  formulation  in  terms  of 
auxiliary  variables  as  outlined  in  Chapter  III  equivalent  to  the  single  boundary  condition 
(G.6)  as: 


where: 


Aj-1  A? 

< h  =  ^ 


(G.7) 


0J  =  0 


This  set  of  conditions  involves  only  first-order  derivatives.  However,  due  to  the  appear¬ 
ance  of  the  ^/-derivative  in  (G.7),  one  cannot  discretize  the  (p3  on  the  boundary  alone. 
Therefore  we  shall  manipulate  (G.7)  in  order  to  get  rid  of  the  //-derivative. 

The  function  h  satisfies  the  dispersive,  advective  wave  equation  (11.29)  in  D.  Since 
the  function  <p1  is  obtained  by  applying  the  linear  (constant  coefficient)  operator  {py  +  Adt/ 
to  h,  it  is  clear  that  <i>\  should  also  satisfy  the  same  equation  in  D.  Further,  since  03  is  ob¬ 
tained  by  applying  the  same  linear  operator  j  —  1  times  to  <1)\.  the  functions  o:j  should  satisfy 
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an  equation  like  (G.3).  Using  the  following  identities: 


dyy&j  (  &y 


C\ 


3+ 1 


9t  I  (  dy  + 


C 


3+ 1 


A )  <j)j  + 


3  '  2 

0i+i 


dyt4*j  (ydy(f)j 


and  combining  with  (G.3)  and  (G.7)  allows  us  to  formulate  the  boundary  as: 


(2V  V2-cl\ 

V  cj  q  J 


h- 1  + 


2UU 


-  A-i  -  (V2  -  <3)  <t>U+ 

-  21/)  4  -  2I/V0/  -  /Vj-i  = 


C v 2  -  co)  <A+i 


for  j  =  1, . . . ,  J  -  1  (G.8) 


In  (G.8)  and  elsewhere,  a  prime  indicates  differentiation  with  respect  to  x  along  T N,  i.e., 
the  tangential  derivative  along  T  \<.  As  desired,  the  new  boundary  condition  (G.8)  does  not 
involve  ^-derivatives.  In  addition,  there  are  no  high-rc  or  t  derivatives  beyond  second  order. 
The  new  formulation  of  the  Jth- order  NRBC  on  T N  can  be  summarized  as  follows: 


Poh  +  dyh  —  <p\  , 

i  3"  i’j—i  T  (3j4>j  70j  /  i  ^y4>j+ 1  (G.9) 

0o  =  h  0j  =  0 


where 


A)  — 


A-  =  (^2  - 


K'j  = 


2C/U 


-217, 


7  =  2C/V, 


A,  =  U2  -  c2, 
A  X  =  U2-  cl 
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c. 


FORMULATION  FOR  THE  SOUTHERN  BOUNDARY 


The  Higdon  boundary  condition  for  Ts  is  given  by: 


0  on  r5. 


(G.10) 


When  imposed  as  a  boundary  condition  on  T 5,  we  can  recast  this  formulation  in  terms  of 
auxiliary  variables  as  outlined  in  Chapter  III  equivalent  to  the  single  boundary  condition 
(G.  10)  as: 


Cs‘) 


where:  4>0  =  h 


j  =  (G.ll) 


<t>J  =  0 


This  set  of  conditions  involves  only  first-order  derivatives.  However,  due  to  the  appear¬ 
ance  of  the  ^-derivative  in  (G.ll),  one  cannot  discretize  the  0j  on  the  boundary  T^  alone. 
Therefore  we  shall  manipulate  (G.l  1)  in  order  to  get  rid  of  the  ^-derivative. 

The  function  h  satisfies  the  dispersive,  advective  wave  equation  (11.29)  in  I).  Since 
the  function  cf>i  is  obtained  by  applying  the  linear  (constant  coefficient)  operator  (dy  ~  Ad\ 
to  h,  it  is  clear  that  0 1  should  also  satisfy  the  same  equation  in  I).  Further,  since  <fij  is  ob¬ 
tained  by  applying  the  same  linear  operator  j  —  1  times  to  <t)\.  the  functions  o:)  should  satisfy 
an  equation  like  (G.3).  Using  the  following  identities: 

dyy(f)j 

dyt^j 
^xy 


°,  +  c 


-a 


3+ 1 


d„  —  — — <9^  (j)j  + 


C 


3+ 1 


r  2  ^ 
uj+ 1 


d t  (dyfa) 
dx  (dyCpj) 
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and  combining  with  (G.3)  and  (G.ll)  allows  us  to  formulate  the  boundary  as: 


(  w  1 

V  Cj  c]  ) 


7-i  -  +  2 uj  ^  -  (t/2  -  c2)  4>U- 

(  A  +  ch)  +  2V)  *  "  2UV4li  ~  /2*"‘  =  -  c»)  h*1 


for  j  =  (G.12) 


In  (G.12)  and  elsewhere,  a  prime  indicates  differentiation  with  respect  to  x  along  i.e., 
the  tangential  derivative  along  IV  As  desired,  the  new  boundary  condition  (G.12)  does  not 
involve  ^-derivatives.  In  addition,  there  are  no  high- a;  or  t  derivatives  beyond  second  order. 
The  new  formulation  of  the  .V'-ordcr  NRBC  onfs  can  be  summarized  as  follows: 


l%h  -  dyh  =  , 

1  +  —  A xftj-l  +  PjPj  —  7 ~~  /20j- 1  =  Ay0j  +  1 


J 

=  h 


(G.13) 


0j  =  0 


where 


/ %  — 
Pj  = 


1  2V 

"cp 


V 2 


q 


Kj  = 


2  UV 

~c~ 


2 U,  Xy  =  V2  -  cl 


72  ~  (7  +  77)  “ 2V’  7  =  2UV’ 


K  ~  U  -  4 
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APPENDIX  H.  OPEN  PLANE  DOMAIN  ROTATION  IN  THE 
DIRECTION  OF  ADVECTION 


Here,  we  consider  the  effects  of  diagonal  advection  and  how  our  formulation  may 
be  simplified.  Suppose  we  have  an  open  plane  domain  where  NRBCs  are  specified  on 
all  four  cardinal  boundaries.  The  question  arises  -  since  we  are  dealing  with  an  infinite 
domain,  can  the  problem  be  simplified  by  a  change  in  coordinates?  If  so,  how  would  this 
adjust  the  problem,  and  would  it  make  the  problem  any  easier? 

To  examine  this  question,  recall  our  PDE  in  its  standard  x  —  y  plane: 

(dt  +  Udx  +  Vdyf  h  -  CqV2/i  +  fh  =  0  (H.l) 

Further,  suppose  that  the  advection  velocities  U  and  V  are  non-zero  in  both  directions 
resulting  in  activation  of  all  components  of  (H.l).  To  fix  some  ideas,  let  U  >0  and  V  >  0. 
The  goal  is  to  convert  (H.l)  in  its  current  coordinate  system  to  one  that  is  in  the  direction 
of  the  advection  velocity.  To  see  this,  consider  Figure  35. 


7/  y 


Figure  35:  Generation  of  a  new  coordinate  system  in  the  direction  of  advection 
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The  new  coordinate  system  (£,  rj)  places  the  the  £  axis  in  the  direction  of  advection. 
This  transformation  is  a  simple  rotation  that  can  be  described  by: 

x  =  cos(0)£  —  sm(9)rj 
y  =  sin(#)£  +  cos(9)r] 


or: 


£  =  cos  {9)x  +  sin  {6)y 
y  =  —  sin(6))a;  +  cos  {9)y 


Clearly  U  and  V  are  also  related  by  the  geometry  of  the  problem. 


sin(0) 

cos(9) 


V 

Vu 2  +  c2 
u 

Vu2  +  c2 


(H.2) 


We  note  now  that  we  can  express  h  in  terms  of  the  new  coordinate  system  as  h(x,  y,  t)  = 
h(£(x,  y),r)(x,  y),t ).  Since  (H.l)  is  developed  in  the  (x,  y)  system,  we  use  the  chain  rule 
to  expand  (H.l)  in  terms  of  (£,  rj).  We  adopt  the  shorthand  convention 


_  dli  _  d2h 

''a  ^ ab  77T 

oa  oaob 
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to  yield  the  following  expansions 


h  y  "4  Vy 

hxy  ^l&^xCy  4"  h^i i  i^r  fjy  ~f  ^yVx )  hrjy  ^jx^jy  H-  h^C,Xy  3“  hyTjxy  (H.3) 

hxx  4^4  2  h^r/^xYjx  4"  h"ipiVx  4"  h^Cxx  hjjTjx x 

h yy  +  2h^yTjy  -\-  hvvr)y  “h  44?/  4"  hy^Jyy 

Since  the  coordinate  transformation  is  linear,  the  problem  is  simplified  even  more 
as  the  second  order  metric  terms  vanish. 

=  cos  0  rjx  —  —  sin  0  4  =  sin  6  rjy  —  cos(9) 

(H.4) 

4a:  4.V  4y  9xx  9yy  9xy  0 

Now,  looking  term  by  term  at  (H.l),  in  light  of  (H.3),  we  consolidate  terms  to  see 

the  result  of  the  transformation. 

h  +  Ah#  +  B  +  C  h,£t  +  D  hvt  +  E  h^v  +  f2h  =  0  (H.5) 

where  h  =  44  77,  t)  and: 

A=(U2-  cl)  g  +  (V2  -  c20)  g  +  2 UVtey 
M=(U2-  cl)  r?x  +  (V2  -  eg)  r?2  +  2C/FW 

C  =  2  (UCx  +  VQ  (H.6) 

D  =  2  ( Urjx  +  Vrjy) 

e  =  (f/2  -  cl)  2ixVx  +  (y2  ~  Co)  24^  +  2£/y  (4^  +  4^) 

We  continue  by  using  the  information  provided  by  the  geometry  of  the  transforma¬ 
tion  in  (H.2)  and  the  metric  terms  found  in  (H.4).  Specifically,  if  we  define  the  adjusted 
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advection  velocities  U  and  V  corresponding  to  the  £  and  //  directions  respectively  as: 


u  =  uzx  + 1/4 

V"  =  Ur}x  +  Vrjy 

then  each  of  the  terms  in  (H.6)  simplify  to: 

A  =  U2-  Cl 
B  =  V2  -  Cl 

C  =  2U  (H.7) 

B  =  2V 

E  =  2UV  -  2cl(2£xr)x  +  £ yrjy ) 

’v‘  ^ 

=0 


Of  course,  if  one  examines  V  we  find: 


V  =  Ur)x  +  Vrjy  =  —U  sin(0)  +  V  cos{6)  =  -U 


V 


+  v^- 


u 


Vu 2  +  v2  vc2  +  v2 


=  0  (H.8) 


This  significantly  simplifies  the  problem  to: 


htt  +  {U~  —  Cl )  h ££  —  clhvv  +  2 Uh^t  +  f2h  —  0 


(H.9) 


This  procedure  shows  that  for  the  open  plane  domain,  any  constant  advection  ve¬ 
locity  not  in  a  cardinal  direction  in  the  standard  x  —  y  coordinate  system  can  be  converted  to 
an  equivalent  problem  where  the  advection  velocity  is  in  a  cardinal  direction  in  an  alternate 
£  —  rj  coordinate  system.  The  result  is  that  when  examining  the  open  plane  domain,  one 
only  needs  to  examine  advections  in  one  cardinal  direction  since  a  problem  with  diagonal 
advection  could  be  recast  into  a  cardinal  direction  advection  problem  in  another  coordinate 
system. 


132 


APPENDIX  I.  ARBITRARILY  SHAPED  BOUNDARY 
FORMULATION 


Recall  the  KGE 


{dt  +  Udx  +  Vdyf  h  -  c20V2h  +  f2h  =  0 


(I-D 


for  which  we  wish  to  formulate  the  G-N  auxiliary  variable  formulation  for  an  arbitrarily 
shaped  boundary.  The  Higdon  boundary  condition  of  order  J  is  given  by 


Hj: 


h  =  0 


on  r 


(1.2) 


Now,  we  introduce  the  auxiliary  functions  0 i, Oj-\,  which  are  defined  on  T  as  well  as 
in  the  exterior  domain  D  (see  Figure  29).  Eventually,  we  shall  use  these  functions  only  on 
T,  but  the  derivation  requires  that  they  be  defined  in  I)  as  well,  or  at  least  in  a  non- vanishing 
region  adjacent  to  T.  The  functions  0j  are  defined  via  the  relations 

(dn  +  h  =  ,  (1.3) 

01  =  02  r  (1-4) 


(dn  +  0j_i  =  0  .  (1.5) 

By  definition,  these  relations  hold  in  D,  and  also  on  T.  It  is  easy  to  see  that  (1.3  - 1.5),  when 
imposed  as  boundary  conditions  on  T,  are  equivalent  to  the  single  boundary  condition  (1.2). 
If  we  also  define 

0o  =  h  0j  =  0  ,  (1.6) 
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then  we  can  write  (1.3  - 1.5)  concisely  as 


(dn  +  (j)j- 1  =  (pj  j  =  1, . . . ,  J  .  (1.7) 

This  set  of  conditions  involves  only  first-order  derivatives.  However,  due  to  the  appearance 
of  the  normal  derivative  in  (1.7),  one  cannot  discretize  the  <pj  on  the  boundary  T  alone. 
Therefore,  we  shall  manipulate  (1.7)  in  order  to  get  rid  of  the  normal  derivative  dn  = 
nxdx  +  riydy  where  nx  =  nx(x,  y )  and  ny  =  ny(x,  y ). 

The  function  h  satisfies  the  KGE  in  I).  Since  the  function  (p1  is  obtained  by  applying 
the  linear  operator  [dn  +  ^rdt^j  to  h,  we  must  consider  what  equation  that  0\  satisfies.  We 
begin  by  applying  this  operator  to  the  KGE  to  yield 

i^dn  +  (j1  ^3'dXxJl  +  A ydyyh  +  2 U dxth  +  2V dyth  +  2 UV dxyh  +  f2h^j  =  0 

^ x  {^xtt^  +  ^■x&xxxh  T  ^ydXyyh  2Udxx^h  2Vdxyth  -t-  A^O^y/i  f  4)xh^ 

Tty  ^ dyt.t.h  ^x4txxyh  Ay<9yyy(^  2 U Dxyth  T  2V dyyth  \xOXyyh  f  OyhPj  (E8) 

7T  ( dttt.h  +  A x.dxxt.h  +  A ydyyth  +  2Udxtth  +  2 Vdytth  +  A xdxyth  +  f  2dthP)  =  0 


Now,  consider  the  first  order  derivatives  for  the  0\  boundary  condition 


9x4^i 

dy(p  1 


A 

Ox 

Ft 


( nx )  4)xVj  T  fix^xxh  ~t~ 
4)xu  T  Tixdxyh  T 
(n*)  0xu  T  nxdxth  -|- 


^  (^y)  4)yU  T  nydxyh  +  ^  r  dxfh 
d  1 

—  (ny)  (9yW  +  Tlydyyh  +  —^dyth 

d  1 

. . ;  (;;?/ )  dyii  T  fi,j Oy i  h  ~\~  . ,  (hi  h 
Ot  L  i 


Since  the  components  of  the  normal  vector  are  themselves  functions  of  a;  and  y,  we  incur 
additional  terms  for  the  derivatives  of  those  components.  If  the  boundary  is  fixed  in  time, 
clearly  the  time  derivatives  of  the  normal  components  will  vanish,  however,  the  spatial 
derivatives  will  not.  We  now  consider  a  simplifying  assumption  that  the  curvature  on  the 
boundary  is  negligible.  In  other  words,  the  spatial  rate  of  change  of  the  components  of  the 
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normal  vector  are  very  small  in  comparison  to  the  other  terms.  This  allows  us  to  neglect 
the  following  terms: 


dxM  ~  dx^y)  ~  dyM  ~  dy  (ny)  “  °’ 


Using  this  assumption,  the  second  order  derivatives  of  fa  are 


&xx4*l  Flxdxxxh  T  Ylydxxyh  "T  s  r  ^xxt^ 

t'l 

&xy 0 1  U/: 9XXy ^  ^y9xyyh  T“  .  ,  Cry/ ^ 

Or 

dyytftl  r/j: Cryy h  T*  ^y9yyy ^  T  .  ,  dyyth 

t  -  1 

9tt.*Pi  Oxtt h  T ■  fiydytth  +  . ,  C// h 

1 

9xt0 i  9XXfh  T  iT'ydxyth  T  . ,  dxtt h 

v  1 

9yt.(p  1  /r.y: Oryt h  +  Tlydyyth  T  ,  ,  9ytt ^ 

L-i 


(1.9) 


We  can  now  use  (1.9)  in  (1.8)  to  find  that  0!  should  also  satisfy  the  KGE  in  I).  Further, 
since  (f)j  is  obtained  by  applying  the  same  operator  to  0j_i,  the  functions  0-j  should  satisfy 
a  similar  equation,  namely, 


<Ct  +  ( U 2  —  Cq)  +  (U2  —  Co)  dyy  +  2C7 <9^  +  2U <9yt  +  2UV dxy  +  f~ 


=  0  (1.10) 


We  now  note  that  the  boundary  condition  and  the  PDE  that  cj)j  satisfies  are  described 
in  two  different  coordinate  systems  -  namely  (n,  r)  and  (x,  y)  respectively.  Consider  an 
arbitrary  part  of  the  boundary  (T)  shown  in  Figure  30.  Of  course,  in  the  most  general 
case,  the  normal  and  tangential  vectors  are  dependent  on  the  position  on  the  boundary,  but 
can  be  computed  given  a  particular  domain.  Since  these  components  then  are  “known,” 
we  consider  a  change  of  coordinates  from  x  and  y  to  n  and  r  as  defined  by  the  linear 
transformation  and  its  associated  differentiation  operator  (VII. 2).  Rewriting  each  operator 
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in  (1. 10)  and  again  using  the  small  curvature  assumption  yields: 


dtt  =  dtt 

dXX  H'xdnn  —  lly  1 1  y  ()r/  r  +  nydTT 

dyy  ttydnn  2tIx'TI 'ydnr  “t“ 

d'd  Ylxdnt  Tlydq-t 

dyt  ^  yd„t  +  nxdTt 

dxy  ^’X^ydnn  ~f  x  )  dnT  TlxTlydTT • 

Substituting  these  terms  back  into  (1. 10)  and  organizing  this  simplifies  to: 

[dtt  +  Adnn  +  BdTT  +  Cdl  +  2  Bdnt  +  2  EdTt  +  f2]  fa  =  0  (1.11) 

where: 

A=(U2-  c20)  n2x  +  ( V 2  -  Cl)  n2y  +  2 UVnxny 
B  =  (U2  —  c20)  n2y  +  (V2  -  Cl)  n2x  -  2 UVnxny 
C  =  2  ( -U 2  +  V2)  nxny  +  2UV  (n2x  -  n2y) 

D  =  Unx  +  Vny 
E  =  — Uny  +  Vnx 

We  now  see  that  this  formulation  is  expressed  in  terms  of  the  normal  and  tangential  coordi¬ 
nate  system  and  further  has  the  same  form  as  (III.  17).  We  can  proceed  in  the  same  manner 
as  outlined  in  Chapter  III  to  eliminate  the  normal  derivatives  to  yield  the  new  formulation 
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of  the  Jth- order  NRBC  on  F  summarized  as  follows: 


M + §; 

<*3<i>j- 1  +  Ij&j-i  ~  B  0?- 1  +  Pifa  ~  c  0|  ~  /20j-i 


^"“cy 


li  -  7S  2  E>  Pi  —  A  7T  + 


j 

2D  A 

“  ~ci  ~  1  ~~  67?’ 


j 


Cf  Ci 


-2D 


3  ^3+ 1  / 

</>o  =  rj  (j)j  =  0. 


01  ) 

A  0j+1 


(1.12) 

(1.13) 


(1.14) 


Here,  prime  indicates  tangential  differentiation  along  T.  As  desired,  the  new  boundary 
condition  does  not  involve  any  normal  derivatives  and  there  are  no  high  tangential  or  time 
derivatives  beyond  second  order. 

There  are  a  few  remaining  concerns.  First,  we  must  be  able  to  compute  the  normal 
and  tangential  components  of  the  vectors  mandated  by  the  mapping  (VII.  2).  Addition¬ 
ally,  we  have  several  terms  which  require  the  integration  of  tangential  derivatives  along  a 
general  boundary.  The  question  arises,  how  do  we  evaluate  these  terms  along  a  particular 
boundary?  Finally,  we  must  still  relate  this  boundary  formulation  back  into  the  interior 
formulation  and  consider  the  appropriate  values  for  the  Cj  terms. 

A.  COMPUTING  THE  NORMAL  AND  TANGENTIAL  VECTOR  COMPONENTS 

Consider  (VII. 2),  which  gives  us  a  way  to  write  the  normal  and  tangential  deriva¬ 
tives  in  terms  of  the  standard  x  —  y  coordinates.  We  can  extend  this  via  the  chain  rule  to 
map  each  of  these  components  in  terms  of  our  canonical  £  —  ;/  coordinates.  Before  we  do 
any  of  this,  however,  we  must  first  be  able  to  compute  the  normal  vectors  at  each  point 
on  the  boundary.  To  see  this,  consider  the  normal  and  tangential  vectors  of  a  canonical 
element  as  shown  in  Figure  36. 
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Figure  36:  Normal  and  tangential  vectors  along  a  canonical  element  boundary. 

After  introducing  the  nonsingular  mapping  £  =  £(x,  y)  and  rj  =  r](x,y),  these 
normal  vectors  [47]  are  given  by: 


n  =  r/ 


|Vj,| 


j?=±1 


n  =  £ 


V£ 


|V£| 


c=±i 


Each  of  these  normal  vectors  can  then  be  used  to  compute  the  tangential  vectors  by 
taking  the  cross  product  of  the  normal  vectors  with  the  unit  vector  (0, 0,  —  1)T.  In  the  case 
of  Figure  36,  the  normalized  tangential  vectors  are: 


When  put  in  terms  of  (IV.8),  the  unit  tangential  vectors  are: 


T\  = 


T3  = 


(1.15) 
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where: 


We  now  have  a  way  to  compute  the  normal  and  tangential  components  required  for  this 
formulation. 

B.  INTEGRATION  OF  TANGENTIAL  DERIVATIVES 


In  (1.13),  we  have  several  terms  that  will  require  the  integration  of  first  and  second 
order  tangential  derivatives  along  a  general  boundary  (after  weak  integral  formulation). 
The  question  arises,  how  do  we  evaluate  these  terms  along  a  particular  boundary? 


1.  Integration  of  First  Order  Tangential  Derivatives 

To  examine  this  arbitrary  domain  formulation  in  greater  detail,  consider  a  single 
element  where  we  evaluate  an  integral  that  contains  a  first  order  tangential  derivative  along 
a  single  canonical  side. 


dT  =  r-V 


Q  =  number  of  quadrature  points  on  a  side 


Here,  J*|  is  the  Jacobian  of  the  transformation  along  the  side  (side  length)  and  wq  is  the 
quadrature  weight  for  Lobatto  integration.  If  we  use  collocated  quadrature,  the  cardinality 
property  of  the  basis  functions  ensures  that  only  the  ith  basis  function  is  nonzero  (in  fact 
unity)  at  quadrature  point  i.  This  rids  us  of  the  summation  over  all  quadrature  points  for 
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basis  function  i  yielding: 


5^)  =w,r,(T,W+TiW 

■'  dr  *  “ciy 


Now,  expand  our  variable  0  in  terms  of  our  canonical  coordinates  (£,  rj)  and  perform  basis 
function  expansions  using  our  2-D  basis  functions  (Mjv  =  (N  +  l)2). 


Mjv 


Mjv 


i=i 

_  ( <9-0*  90  +  <90*  dy 


.7=1  ^ 


*  /  y 

J=i 


9£  <9a:  dq  dx 


Mjv 


j=l 


90]  <90  90]  drf 


9£  9?/  dr]  dy 


Of  course  our  boundary  integral  implies  that  we  are  integrating  strictly  on  a  side.  For 
definiteness  in  the  example,  consider  side  2  where  £  =  +  1 ,  //  e  [—1,  +1].  Now,  we  expand 
using  (IV.8)  and  (1.15): 


Simplifying  and  using  the  definition  of  the  element  Jacobian  (IV.7),  our  task  now  requires 
evaluating: 


=Wi 


/  dx 1  dy 1  dy1  dx 1  \  90]  \ 

\  dy  9£  9/y  9£  /  ^0  dy  3  J 


Mjv 


=  VJ, 


y^90]( 

^  dn 
.7  =  1 


d'tpi 

Now,  consider  the  term  -0,  which  means  “the  partial  derivative  with  respect  to  // 
of  the  jth  basis  function  evaluated  at  the  ith  quadrature  point”.  We  note  that  the  2-D  basis 
functions  ipj  are  formed  using  tensor  products  of  1-D  basis  functions  z/00)  and  z//(9/  j  such 
that  j  =  1, . . . ,  Mn  where  the  mapping  from  1-D  to  2-D  is  j  =  {(k  —  1)(N  +  1)  +  l  : 
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k,  l  —  1, . . .  N  +  1}.  Therefore: 


Using  this  and  the  cardinality  property  of  the  basis  functions,  we  can  evaluate  the  integral 


of  the  first  order  tangential  derivative  as: 


f  ,  d<f>(x,  y )  7-p  duj 


A  similar  argument  holds  for  other  canonical  sides  being  evaluated  on  the  boundaries. 

2.  Integration  of  Second  Order  Tangential  Derivatives 

Recall,  the  second  order  term  from  the  auxiliary  formulation  of  the  form: 


/  d2</>  ^  f  d  (  d<p\  f  d'lpi  d(j) 

dr‘ = X,  r 5F )  dr‘ -Jr,wfr dT‘ 

d<t>  end  f  8 ^  deb 
f  'dr  start  <9r  <9r  ^  s' 

=0  if  closed  boundary 


Now,  suppose  that  we  expand  each  term,  as  in  the  previous  section. 


/‘  ch/y  d(f>  f  (  d^i  d^i\  f  d<f>  d(p\ 

L  ** dr*  =  X.  (T*  &  +  T”  %  J  Xs  +  T^J <r* 


Q 

q=  1 


■ si  Tqzz±.  +  r<7 
V  s  ^  y  By 


d(f)q  d  (f)q 

Txlkc+Tyl)y 
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Consider  ®  evaluated  on  side  2. 


r  Ty 


dx 


dy 


rq  (  d£q  diff*  drf\  q  f  dtfi  d£q  dip?  drf\ 

x  \  d£  dx  dr)  dx  )  Ty  \  d£  dy  dr)  dy  ) 

1  dx  ( d'dl  (  1  dy  \  d'dl  f  1  <9y  \  \ 

\Jq\dr)  V  d£  \\J*\dr,)  +  dy  \  \Jq\  d£j  ) 

1  dy  f  dip?  f  1  dx\  dip?  /  1  Ox \  \ 

Jq\dy  V  ^  V  \Jq\  drU  +  dy  \\J%\dt) ) 

1  /  dxdy  dy  dx  \  dip? 

\Jq\\  Jq I  V  dy  d£  +  dy  d£ )  dy 
1  dj>i 
I  Jq  I  dy 


Now,  combining  this  result  with  the  expansion  of  (2)  as  found  in  Appendix  I.B,  we  find: 


As  in  Appendix  I.B,  similar  results  can  be  obtained  for  each  of  the  other  canonical  sides. 


C.  RELATING  THE  BOUNDARY  AND  INTERIOR  FORMULATIONS 


The  final  four  terms  of  (IV.3)  are  all  boundary  integrals  where  the  NRBC  must  be 
applied.  To  simplify  our  discussion,  we  note  that  these  integrals  only  apply  to  h  on  T,  thus 
we  will  denote  these  terms  as  h.  Much  like  the  auxiliary  variable  formulation,  we  note  that 
the  NRBC  is  defined  in  terms  of  normal  derivatives  (dn)  while  the  boundary  terms  of  (IV.3) 
are  defined  in  terms  of  standard  Cartesian  derivatives  (dx  and  dy).  Again,  we  consider  a 
linear  transformation  of  the  final  four  terms  as  defined  in  (VII. 2).  This  yields,  for  the  first 
boundary  term: 

dr  =  \x  j  % 


dh  f  dh 

T®,  nx  dF  =  Xx  /  To—  nx 

OX  Jr  OX 


dh 


dh 


Ax  rs  fly 

dn  dr 


nx  dT  (1.16) 
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Now,  we  can  directly  use  (1.3)  to  join  the  two  formulations: 


.  dh  dr) 

I  nx- - ny— 

on  or 


nr  dr 


nxdY  (1.17) 


Now,  we  note  that  the  boundary  integral  can  be  discretized  on  the  boundary  alone  -  in 
terms  of  tangential  derivatives,  time  derivatives,  and  the  auxiliary  variable  (f>\.  Performing 
these  substitutions  in  each  boundary  term  of  (IV.3),  we  get  the  revised  weak  form  of  the 
problem: 


^ ih  dft  -  Xx  , 

/o  C/X  C/X 


^  dh  da  -  a. 


d^i  dh 

dll 


dy  dy 


-UV 


O'k  dh 

1  dVt-UV 


dy  dx 


d^i  dh 

dll 


dx  dy 


dh 


dh 


+2 u  /  yt—dn  +  2v  /  %—dn  +  r  /  %hdn 


'  dx 


dy 


dh 

+A*  /  ^i^inl  dT  -  Xx  I  'P i'h  fJ0n2x  dT  -  Xx  /  ^i—nynx  dr 


(1.18) 


+Xy  J  ^i4>\n2y  dr  -  Xy  J  'k ih  (30n2y  dr  +  Xy  J  ^~nxny  dT 
+UV  J  VifanxTiy  dT  —  UV  J  rih/30nxny  dr  —  UV  J  r,j^n2y  dr 

+UV  f  'Ti(j)1nxny  dV  -  UV  f  \k  Ji  /30nxny  dr  +  UV  f  vk,:^n^.  dr  =  0, 

Jr  Jr  Jr  dr 


which,  as  desired,  evaluates  boundary  integrals  with  data  derived  from  only  the  boundary. 


D.  SELECTION  OF  Cj  TERMS 


If  we  examine  the  auxiliary  variable  formulation  (1.13),  we  see  that  the  selection  of 
appropriate  Cj  values  has  yet  to  be  addressed.  As  has  been  previously  shown,  any  choice 
of  Cj  is  guaranteed  to  reduce  spurious  reflection  as  the  order  of  the  NRBC  ( J)  increases. 
Armed  with  this,  we  choose  convenient  values  for  our  Cj  ’s  that  cause  the  second  order  in 
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time  (ay)  terms  to  vanish  for  all  j.  In  this  case: 


„  2D  1  A 

“ 0  “  a  ~ 1  “  c? 


=>  a  =  D  ±  vD  -  A 


=  Unx  +  Vny  ±  y eg  (n^  +  n2y) 
=  Unx  +  Wiy  ±  Co 


This  implies  a  special  value  of  Cj  for  each  boundary  point,  independent  of  NRBC 
order,  i.e.,  C3  =  C(x.  y)  for  all  j.  Further,  this  implies  that  the  other  coefficients  in  the 
auxiliary  formulation  are  independent  of  NRBC  order: 


C 

C(x,y) 


-2E, 


Pj  =  P 


2A 

C(x,y) 


-2D 


The  Higdon  boundary  condition,  while  general  in  nature,  implicitly  assumes  that  by  the 
time  a  wave  pulse  gets  to  the  artificial  boundary,  it  is  traveling  primarily  as  a  plane  wave 
normal  to  the  boundary.  This  choice  for  C3  can  be  thought  of  as  a  choice  that  accounts  for 
any  advection  and  “corrects”  for  the  geometry  of  the  problem  -  i.e.,  non-normal  impinge¬ 
ment  on  the  boundary. 
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